/* nag_opt_check_2nd_deriv (e04hdc) Example Program. * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void hess(Integer n, double xc[], double fhesl[], double fhesd[], Nag_Comm *comm); static void funct(Integer n, double xc[], double *fc, double gc[], Nag_Comm *comm); #ifdef __cplusplus } #endif int main(void) { #define NMAX 4 Integer exit_status=0, i, j, k, n; NagError fail; Nag_Comm comm; double f, *g=0, *hesd=0, *hesl=0, *x=0; INIT_FAIL(fail); #define X(I) x[(I)-1] #define HESL(I) hesl[(I)-1] #define HESD(I) hesd[(I)-1] #define G(I) g[(I)-1] Vprintf("nag_opt_check_2nd_deriv (e04hdc) Example Program Results\n\n"); /* Set up an arbitrary point at which to check the derivatives */ n = NMAX; if (n>=1) { if ( !( hesd = NAG_ALLOC(n, double)) || !( hesl = NAG_ALLOC(n*(n-1)/2, double)) || !( g = NAG_ALLOC(n, double)) || !( x = NAG_ALLOC(n, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } } else { Vprintf("Invalid n.\n"); exit_status = 1; return exit_status; } X(1) = 1.46; X(2) = -0.82; X(3) = 0.57; X(4) = 1.21; Vprintf("The test point is\n"); for (j = 1; j <= n; ++j) Vprintf("%9.4f", X(j)); Vprintf("\n"); /* Check the 1st derivatives */ /* nag_opt_check_deriv (e04hcc). * Derivative checker for use with nag_opt_bounds_deriv * (e04kbc) */ nag_opt_check_deriv(n, funct, &X(1), &f, &G(1), &comm, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_check_deriv (e04hcc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Check the 2nd derivatives */ /* nag_opt_check_2nd_deriv (e04hdc). * Checks second derivatives of a user-defined function */ nag_opt_check_2nd_deriv(n, funct, hess, &X(1), &G(1), &HESL(1), &HESD(1), &comm, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_check_2nd_deriv (e04hdc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\n2nd derivatives are consistent with 1st derivatives.\n\n"); Vprintf("At the test point, funct gives the function value, %12.4e\n", f); Vprintf("and the 1st derivatives\n"); for (j = 1; j <= n; ++j) Vprintf("%12.3e%s", G(j), j%4?"":"\n"); Vprintf("\nhess gives the lower triangle of the Hessian matrix\n"); Vprintf("%12.3e\n", HESD(1)); k = 1; for (i = 2; i <= n; ++i) { for (j = k; j <= k + i - 2; ++j) Vprintf("%12.3e", HESL(j)); Vprintf("%12.3e\n", HESD(i)); k = k + i - 1; } END: if (hesd) NAG_FREE(hesd); if (hesl) NAG_FREE(hesl); if (g) NAG_FREE(g); if (x) NAG_FREE(x); return exit_status; } static void funct(Integer n, double xc[], double *fc, double gc[], Nag_Comm *comm) { /* Routine to evaluate objective function and its 1st derivatives. */ #define GC(I) gc[(I)-1] #define XC(I) xc[(I)-1] *fc = pow(XC(1)+10.0*XC(2), 2.0) + 5.0*pow(XC(3)-XC(4), 2.0) + pow(XC(2)-2.0*XC(3), 4.0) + 10.0*pow(XC(1)-XC(4), 4.0); GC(1) = 2.0*(XC(1)+10.0*XC(2)) + 40.0*pow(XC(1)-XC(4),3.0); GC(2) = 20.0*(XC(1)+10.0*XC(2)) + 4.0*pow(XC(2)-2.0*XC(3),3.0); GC(3) = 10.0*(XC(3)-XC(4)) - 8.0*pow(XC(2)-2.0*XC(3),3.0); GC(4) = 10.0*(XC(4)-XC(3)) - 40.0*pow(XC(1)-XC(4), 3.0); } static void hess(Integer n, double xc[], double fhesl[], double fhesd[], Nag_Comm *comm) { /* Routine to evaluate 2nd derivatives */ #define FHESD(I) fhesd[(I)-1] #define FHESL(I) fhesl[(I)-1] #define XC(I) xc[(I)-1] FHESD(1) = 2.0 + 120.0*pow(XC(1)-XC(4), 2.0); FHESD(2) = 200.0 + 12.0*pow(XC(2)-2.0*XC(3), 2.0); FHESD(3) = 10.0 + 48.0*pow(XC(2)-2.0*XC(3), 2.0); FHESD(4) = 10.0 + 120.0*pow(XC(1)-XC(4), 2.0); FHESL(1) = 20.0; FHESL(2) = 0.0; FHESL(3) = -24.0*pow(XC(2)-2.0*XC(3), 2.0); FHESL(4) = -120.0*pow(XC(1)-XC(4), 2.0); FHESL(5) = 0.0; FHESL(6) = -10.0; }