Bezier and B - spline surface interpolation

 We are given an array pp  of size (nu+llu)(nv+llv), containing data points to be interpolated. Data points have parameter values at the Greville abscissae.

nu=8;  (*u-degree *)
nv=8; (*v-degree*)
llu=1; (* number of u-domain intervals *)
llv=1; (* number of v-domain intervals *)
cols=nv+llv;
rows=nu+llu;
knotu = Join[Table[0,{i,0,nu}],Table[i/llu,{i,1,llu-1}],Table[1,{i,0,nu}]]//N;
knotv = Join[Table[0,{i,0,nv}],Table[i/llv,{i,1,llv-1}],Table[1,{i,0,nv}]]//N;

(* Greville abscissae grevu, grevv: *)
grevu=Table[Sum[knotu[[j+1]],{j,i,i+nu-1}]/nu,{i,1,nu+llu}]//N;
grevv=Table[Sum[knotv[[j+1]],{j,i,i+nv-1}]/nv,{i,1,nv+llv}]//N;

(* Vandermonde matrices vandu, vandv: *)
vandu=Table[BSplineBasis[{nu,knotu},j,grevu[[i]]],{i,1,nu+llu},{j,0,nu+llu-1 }];
vandv=Transpose[Table[BSplineBasis[{nv,knotv},j,grevv[[i]]],{i,1,nv+llv},{j,0,nv+llv-1 }]];

(* rectangular array pp of data points: *)
(* airplane wing: nu, nv, llu, llv: 5,1,1,1*)
pp=Table[{.7Cos[i]+ 0.3Sin[2i],Sin[1.1i]+Cos[i], j},{i,nu+llu},{j,nv+llv}]//N;
pp=Table[{i,j,i^2/5-i+Sin[j]},{i,nu+llu},{j,nv+llv}]//N;
pp=Table[{Sqrt[i]*Cos[j],Sqrt[i]*Sin[j],.5Sqrt[(j^2+i^2)]},{i,1,nu+llu},{j,1,nv+llv}]//N;

MatrixForm[vandu];
MatrixForm[pp];
MatrixForm[Inverse[Transpose[vandv]]];

(* Resulting contol net: *)
intdd=Inverse[vandu].Transpose[Inverse[Transpose[vandv]].Transpose[pp]];
(* Don't understand all the Transpose stuff -- but that makes it work. *)

intnetplot=Graphics3D[{PointSize[Large],Green,Map[Point,pp],Gray,Line[intdd],Line[Transpose[intdd]]}];
intddplot=Graphics3D[{PointSize[Large],Red,Map[Point,intdd]}];
intsurfplot=Graphics3D[BSplineSurface[intdd,SplineKnots->{knotu,knotv},SplineDegree-> {nu,nv}]];

Show[{intddplot,intsurfplot,intnetplot},Boxed->False]