/* test_nuderiv3dg.c   test  nuderiv3dg  with max 28 points     */
/*                     order 1  minimum independent points 4    */
/*                     order 2  minimum independent points 10   */
/*                     order 3  minimum independent points 20   */
/*                     order 4  minimum independent points 35   */
/*                     going beyond minimum points does not help much  */

#include <stdio.h>
#include <stdlib.h>
#include "nuderiv3dg.h"
#include "inverse.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  max
#define max(x,y) ((x)>(y)?(x):(y))

/* 1st order test case */

static double f1(double x, double y, double z)
{
  return x + 2.0*y + 3.0*z + 4.0;
}
static double f1x(double x, double y, double z)
{
  return 1.0;
}
static double f1y(double x, double y, double z)
{
  return 2.0;
}
static double f1z(double x, double y, double z)
{
  return 3.0;
}

/* 2nd order test case */
static double f2(double x, double y, double z)
{
  return x*x + 2.0*x*y + 3.0*x*z + 4.0*y*y +
     5.0*y*z + 6.0*z*z + 7.0*x   + 8.0*y   +
    9.0*z   + 10.0;
}
static double f2x(double x, double y, double z)
{
  return 2.0*x + 2.0*y + 3.0*z + 7.0;
}
static double f2y(double x, double y, double z)
{
  return 2.0*x + 8.0*y + 5.0*z + 8.0;
}
static double f2z(double x, double y, double z)
{
  return 3.0*x + 5.0*y + 12.0*z + 9.0;
}
static double f2xx(double x, double y, double z)
{
  return 2.0;
}
static double f2xy(double x, double y, double z)
{
  return 2.0;
}
static double f2xz(double x, double y, double z)
{
  return 3.0;
}
static double f2yy(double x, double y, double z)
{
  return 8.0;
}
static double f2yz(double x, double y, double z)
{
  return 5.0;
}
static double f2zz(double x, double y, double z)
{
  return 12.0;
}


/* 3rd order test case */
static double f3(double x, double y, double z)
{
  return x*x*x +  2.0*x*x*y +  3.0*x*x*z + 4.0*x*y*y +
     5.0*x*y*z +  6.0*x*z*z +  7.0*y*y*y + 8.0*y*y*z +
     9.0*y*z*z + 10.0*z*z*z + 11.0*x*x   + 12.0*x*y +
    13.0*x*z   + 14.0*y*y   + 15.0*y*z   + 16.0*z*z  +
    17.0*x     + 18.0*y     + 19.0*z     + 20.0;
}
static double f3x(double x, double y, double z)
{
  return 3.0*x*x +  4.0*x*y +  6.0*x*z +  4.0*y*y +
         5.0*y*z +  6.0*z*z +
                              22.0*x   + 12.0*y   +
        13.0*z  +
    17.0;
}
static double f3y(double x, double y, double z)
{
  return             2.0*x*x  +             8.0*x*y +
         5.0*x*z +              21.0*y*y + 16.0*y*z +
         9.0*z*z +                         12.0*x   +
                    28.0*y    + 15.0*z   +
    18.0;
}
static double f3z(double x, double y, double z)
{
  return                       3.0*x*x +
         5.0*x*y + 12.0*x*z +             8.0*y*y +
        18.0*y*z + 30*z*z   +
        13.0*x   +            15.0*y   + 32.0*z   +
    19.0;
}
static double f3xx(double x, double y, double z)
{
  return 6.0*x + 4.0*y +  6.0*z + 22.0;
}
static double f3xy(double x, double y, double z)
{
  return 4.0*x + 8.0*y + 5.0*z + 12.0  ;
}
static double f3xz(double x, double y, double z)
{
  return 6.0*x + 5.0*y + 12.0*z + 13.0;
}
static double f3yy(double x, double y, double z)
{
  return 8.0*x + 42.0*y + 16.0*z + 28.0;
}
static double f3yz(double x, double y, double z)
{
  return 5.0*x + 16.0*y + 18.0*z + 15.0;
}
static double f3zz(double x, double y, double z)
{
  return 12.0*x + 18.0*y + 60.0*z + 32.0;
}
static double f3xxx(double x, double y, double z)
{
  return 6.0;
}
static double f3xxy(double x, double y, double z)
{
  return 4.0;
}
static double f3xxz(double x, double y, double z)
{
  return 6.0;
}
static double f3xyy(double x, double y, double z)
{
  return 8.0;
}
static double f3xyz(double x, double y, double z)
{
  return 5.0;
}
static double f3xzz(double x, double y, double z)
{
  return 12.0;
}
static double f3yyy(double x, double y, double z)
{
  return 42.0;
}
static double f3yyz(double x, double y, double z)
{
  return 16.0;
}
static double f3yzz(double x, double y, double z)
{
  return 18.0;
}
static double f3zzz(double x, double y, double z)
{
  return 60.0;
}

/* 4th order test case  TBD*/


int main(int argc, char *argv[])
{
  int i, j;
  int debug = 0;
  int n, npoint;
  int order;
  double x[35] = {0.05, 0.14, 0.27, 0.42, 0.57, 0.71, 0.86,
                  0.04, 0.18, 0.31, 0.43, 0.63, 0.82, 0.95,
                  0.07, 0.16, 0.33, 0.41, 0.59, 0.82, 0.99,
                  0.11, 0.21, 0.35, 0.45, 0.65, 0.75, 0.97,
                  0.06, 0.17, 0.32, 0.46, 0.61, 0.79, 0.98};
  double y[35] = {0.06, 0.05, 0.07, 0.10, 0.13, 0.15, 0.17,
                  0.30, 0.33, 0.36, 0.39, 0.31, 0.32, 0.37,
                  0.48, 0.49, 0.52, 0.51, 0.47, 0.50, 0.53,
                  0.65, 0.77, 0.71, 0.67, 0.67, 0.75, 0.76,
                  0.92, 0.97, 0.93, 0.98, 0.87, 0.91, 0.99};
  double z[35] = {0.07, 0.28, 0.47, 0.69, 0.91,
                  0.15, 0.45, 0.68, 0.93, 0.05,
                  0.26, 0.50, 0.97, 0.08, 0.14,
                  0.39, 0.65, 0.76, 0.93, 0.27,
                  0.56, 0.82, 0.99, 0.11, 0.33,
                  0.72, 0.94, 0.09, 0.22, 0.65,
                  0.95, 0.10, 0.24, 0.42, 0.82};
  double U[35];
  int ip[35];
  double est, actual, dx, dy, dz;
  double cx[35], cy[35], cz[35], cxx[35], cxy[35], cxz[35];
  double cyy[35], cyz[35], czz[35], cxxx[35], cxxy[35];
  double cxxz[35], cxyy[35], cxyz[35], cxzz[35], cyyy[35], cyyz[35];
  double cyzz[35], czzz[35];
  double error, maxerr;
  int point;
  int m = 0;

  printf("test_nuderiv3dg.c running \n");
  self_test(); /* nuderiv3dg.c */
  n = 35;
  npoint = 6;
  for(i=0; i<npoint; i++)
  {
    U[i] = f1(x[i],y[i],z[i]);
    dx = f1x(x[i],y[i],z[i]);
    dy = f1y(x[i],y[i],z[i]);
    dz = f1z(x[i],y[i],z[i]);
    if(debug)
      printf("point %2d  x=%5.3f, y=%5.3f, z=%5.3f, U=%7.3f, dx=%f, dy=%f, dz=%f \n",
	     i, x[i], y[i], z[i], U[i], dx, dy, dz);
  } 

  /* test f1 dx, dy dz */
  order = 1;
  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f1 dx */
  {
    m = nu3dx(order, npoint, point, x, y, z, cx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
    actual = f1x(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f1 dx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f1 dx */
  printf("f1 dx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f1 dy */
  {
    m = nu3dy(order, npoint, point, x, y, z, cy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
    actual = f1y(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug)printf("f1 dy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f1 dy */
  printf("f1 dy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f1 dz */
  {
    m = nu3dz(order, npoint, point, x, y, z, cz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cz[i]*U[ip[i]];
    actual = f1z(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f1 dz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f1 dz */
  printf("f1 dz max error=%g \n", maxerr);

  /* test f2 dx, dy, dz, dxx, dxy, dxz, dyy, dyz, dzz */
  npoint = 12;
  for(i=0; i<npoint; i++)  U[i] = f2(x[i],y[i],z[i]);
  order = 2;
  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dx */
  {
    m = nu3dx(order, npoint, point, x, y, z, cx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
    actual = f2x(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dx */
  printf("f2 dx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dy */
  {
    m = nu3dy(order, npoint, point, x, y, z, cy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
    actual = f2y(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f2 dy */
  printf("f2 dy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dz */
  {
    m = nu3dz(order, npoint, point, x, y, z, cz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cz[i]*U[ip[i]];
    actual = f2z(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f2 dz */
  printf("f2 dz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dxx */
  {
    m = nu3dxx(order, npoint, point, x, y, z, cxx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxx[i]*U[ip[i]];
    actual = f2xx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dxx */
  printf("f2 dxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dxy */
  {
    m = nu3dxy(order, npoint, point, x, y, z, cxy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxy[i]*U[ip[i]];
    actual = f2xy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dxy */
  printf("f2 dxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dxz */
  {
    m = nu3dxz(order, npoint, point, x, y, z, cxz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxz[i]*U[ip[i]];
    actual = f2xz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dxz */
  printf("f2 dxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dyy */
  {
    m = nu3dyy(order, npoint, point, x, y, z, cyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyy[i]*U[ip[i]];
    actual = f2yy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dyy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f2 dyy */
  printf("f2 dyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dyz */
  {
    m = nu3dyz(order, npoint, point, x, y, z, cyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyz[i]*U[ip[i]];
    actual = f2yz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dyz %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f2 dyz */
  printf("f2 dyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dzz */
  {
    m = nu3dzz(order, npoint, point, x, y, z, czz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += czz[i]*U[ip[i]];
    actual = f2zz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f2 dzz */
  printf("f2 dzz max error=%g \n", maxerr);


  /* test f3 dx, dy, dz, dxx, dxy, dxz, dyy, dyz, dzz */
  /*         dxxx, dxxy, dxxz, dxyy, dxyz, dxzz, dyyy, dyyz, dyzz, dzzz */
  npoint = 22;
  for(i=0; i<npoint; i++)  U[i] = f3(x[i],y[i],z[i]);
  order = 3;
  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dx */
  {
    m = nu3dx(order, npoint, point, x, y, z, cx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
    actual = f3x(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dx */
  printf("f3 dx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dy */
  {
    m = nu3dy(order, npoint, point, x, y, z, cy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
    actual = f3y(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dy */
  printf("f3 dy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dz */
  {
    m = nu3dz(order, npoint, point, x, y, z, cz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cz[i]*U[ip[i]];
    actual = f3z(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dz */
  printf("f3 dz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxx */
  {
    m = nu3dxx(order, npoint, point, x, y, z, cxx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxx[i]*U[ip[i]];
    actual = f3xx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxx */
  printf("f3 dxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxy */
  {
    m = nu3dxy(order, npoint, point, x, y, z, cxy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxy[i]*U[ip[i]];
    actual = f3xy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxy */
  printf("f3 dxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxz */
  {
    m = nu3dxz(order, npoint, point, x, y, z, cxz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxz[i]*U[ip[i]];
    actual = f3xz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxz */
  printf("f3 dxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyy */
  {
    m = nu3dyy(order, npoint, point, x, y, z, cyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyy[i]*U[ip[i]];
    actual = f3yy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyy */
  printf("f3 dyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyz */
  {
    m = nu3dyz(order, npoint, point, x, y, z, cyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyz[i]*U[ip[i]];
    actual = f3yz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyz %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyz */
  printf("f3 dyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dzz */
  {
    m = nu3dzz(order, npoint, point, x, y, z, czz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += czz[i]*U[ip[i]];
    actual = f3zz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dzz */
  printf("f3 dzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxxx */
  {
    m = nu3dxxx(order, npoint, point, x, y, z, cxxx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
    actual = f3xxx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxxx */
  printf("f3 dxxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxxy */
  {
    m = nu3dxxy(order, npoint, point, x, y, z, cxxy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
    actual = f3xxy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxxy */
  printf("f3 dxxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxxz */
  {
    m = nu3dxxz(order, npoint, point, x, y, z, cxxz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxz[i]*U[ip[i]];
    actual = f3xxz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxxz */
  printf("f3 dxxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxyy */
  {
    m = nu3dxyy(order, npoint, point, x, y, z, cxyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
    actual = f3xyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxyy */
  printf("f3 dxyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxyz */
  {
    m = nu3dxyz(order, npoint, point, x, y, z, cxyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyz[i]*U[ip[i]];
    actual = f3xyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxyz */
  printf("f3 dxyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxzz */
  {
    m = nu3dxzz(order, npoint, point, x, y, z, cxzz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxzz[i]*U[ip[i]];
    actual = f3xzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dxzz */
  printf("f3 dxzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyyy */
  {
    m = nu3dyyy(order, npoint, point, x, y, z, cyyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
    actual = f3yyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyyy */
  printf("f3 dyyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyyz */
  {
    m = nu3dyyz(order, npoint, point, x, y, z, cyyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyz[i]*U[ip[i]];
    actual = f3yyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyyz %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyyz */
  printf("f3 dyyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dzzz */
  {
    m = nu3dzzz(order, npoint, point, x, y, z, czzz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += czzz[i]*U[ip[i]];
    actual = f3zzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dzzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dzzz */
  printf("f3 dzzz max error=%g \n", maxerr);


  printf("end test_nuderiv3dg.c\n");
  return 0;
} /* end  test_nuderiv3dg.c */