/* 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 #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL h(Integer n, double xc[], double fhesl[], double fhesd[], Nag_Comm *comm); static void NAG_CALL funct(Integer n, double xc[], double *fc, double gc[], Nag_Comm *comm); #ifdef __cplusplus } #endif int main(int argc, char *argv[]) { FILE *fpout; 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); /* Check for command-line IO options */ fpout = nag_example_file_io(argc, argv, "-results", NULL); #define X(I) x[(I) -1] #define HESL(I) hesl[(I) -1] #define HESD(I) hesd[(I) -1] #define G(I) g[(I) -1] fprintf(fpout, "nag_opt_check_2nd_deriv (e04hdc) Example Program Results\n\n"); /* Set up an arbitrary point at which to check the derivatives */ n = 4; 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))) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } } else { fprintf(fpout, "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; fprintf(fpout, "The test point is\n"); for (j = 1; j <= n; ++j) fprintf(fpout, "%9.4f", X(j)); fprintf(fpout, "\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) { fprintf(fpout, "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, h, &X(1), &G(1), &HESL(1), &HESD(1), &comm, &fail); if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_opt_check_2nd_deriv (e04hdc).\n%s\n", fail.message); exit_status = 1; goto END; } fprintf(fpout, "\n2nd derivatives are consistent with 1st derivatives.\n\n"); fprintf(fpout, "At the test point, funct gives the function value, %13.4e\n", f); fprintf(fpout, "and the 1st derivatives\n"); for (j = 1; j <= n; ++j) fprintf(fpout, "%12.3e%s", G(j), j%4?"":"\n"); fprintf(fpout, "\nh gives the lower triangle of the Hessian matrix\n"); fprintf(fpout, "%12.3e\n", HESD(1)); k = 1; for (i = 2; i <= n; ++i) { for (j = k; j <= k + i - 2; ++j) fprintf(fpout, "%12.3e", HESL(j)); fprintf(fpout, "%12.3e\n", HESD(i)); k = k + i - 1; } END: if (fpout != stdout) fclose(fpout); 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 NAG_CALL funct(Integer n, double xc[], double *fc, double gc[], Nag_Comm *comm) { /* Routine to evaluate objective function and its 1st derivatives. */ *fc = pow(xc[0]+10.0*xc[1], 2.0) + 5.0*pow(xc[2]-xc[3], 2.0) + pow(xc[1]-2.0*xc[2], 4.0) + 10.0*pow(xc[0]-xc[3], 4.0); gc[0] = 2.0*(xc[0]+10.0*xc[1]) + 40.0*pow(xc[0]-xc[3], 3.0); gc[1] = 20.0*(xc[0]+10.0*xc[1]) + 4.0*pow(xc[1]-2.0*xc[2], 3.0); gc[2] = 10.0*(xc[2]-xc[3]) - 8.0*pow(xc[1]-2.0*xc[2], 3.0); gc[3] = 10.0*(xc[3]-xc[2]) - 40.0*pow(xc[0]-xc[3], 3.0); } static void NAG_CALL h(Integer n, double xc[], double fhesl[], double fhesd[], Nag_Comm *comm) { /* Routine to evaluate 2nd derivatives */ fhesd[0] = 2.0 + 120.0*pow(xc[0]-xc[3], 2.0); fhesd[1] = 200.0 + 12.0*pow(xc[1]-2.0*xc[2], 2.0); fhesd[2] = 10.0 + 48.0*pow(xc[1]-2.0*xc[2], 2.0); fhesd[3] = 10.0 + 120.0*pow(xc[0]-xc[3], 2.0); fhesl[0] = 20.0; fhesl[1] = 0.0; fhesl[2] = -24.0*pow(xc[1]-2.0*xc[2], 2.0); fhesl[3] = -120.0*pow(xc[0]-xc[3], 2.0); fhesl[4] = 0.0; fhesl[5] = -10.0; }