begin
  NX = 21
  NY = 21
  NZ = 21
  x = (/0.00, 1.00, 0.00, 1.00, 0.40, 0.75/)
  y = (/0.00, 0.00, 1.00, 1.00, 0.20, 0.65/)
  z = (/0.00, 0.00, 0.00, 0.00, 1.25, 0.80/)
  u = x*x + y*y + z*z

  xo = new((/NX/),float)
  yo = new((/NY/),float)
  zo = new((/NZ/),float)

;
; Create the output grid.
;
  xmin = -2.0
  ymin = -2.0
  zmin = -2.0
  xmax =  2.0
  ymax =  2.0
  zmax =  2.0

  ii = fspan(0,20.,21)
  xo = xmin + (ii/(NX-1)) * (xmax-xmin)
  yo = ymin + (ii/(NY-1)) * (ymax-ymin)
  zo = zmin + (ii/(NZ-1)) * (zmax-zmin)

;
; Interpolate.  "output" will be a 3-dimensional grid dimensioned 21x21x21.
;
  output = dsgrid3s(x, y, z, u, xo, yo, zo)
  print(output)
end