/* nag_opt_nlin_lsq (e04unc) Example Program. * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * * Mark 6 revised, 2000. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void objfun(Integer m, Integer n, double x[], double f[], double fjac[], Integer tdfjac, Nag_Comm *comm); static void confun(Integer n, Integer ncnlin, Integer needc[], double x[], double conf[], double cjac[], Nag_Comm *comm); static void user_print(const Nag_Search_State *st, Nag_Comm *comm); #ifdef __cplusplus } #endif static int ex1(void); static int ex2(void); static void objfun(Integer m, Integer n, double x[], double f[], double fjac[], Integer tdfjac, Nag_Comm *comm) { #define FJAC(I,J) fjac[(I)*tdfjac + (J)] /* Initialized data */ static double a[44] = { 8.0, 8.0, 10.0, 10.0, 10.0, 10.0, 12.0, 12.0, 12.0, 12.0, 14.0, 14.0, 14.0, 16.0, 16.0, 16.0, 18.0, 18.0, 20.0, 20.0, 20.0, 22.0, 22.0, 22.0, 24.0, 24.0, 24.0, 26.0, 26.0, 26.0, 28.0, 28.0, 30.0, 30.0, 30.0, 32.0, 32.0, 34.0, 36.0, 36.0, 38.0, 38.0, 40.0, 42.0 }; /* Local variables */ double temp; Integer i; double x0, x1, ai; /* Function to evaluate the objective subfunctions * and their 1st derivatives. */ x0 = x[0]; x1 = x[1]; for (i = 0; i < m; ++i) { ai = a[i]; temp = exp(-x1 * (ai - 8.0)); /* Evaluate objective subfunction f(i+1) only if required */ if (comm->needf == i+1 || comm->needf == 0) f[i] = x0 + (.49 - x0) * temp; if (comm->flag == 2) { FJAC(i,0) = 1.0 - temp; FJAC(i,1) = -(.49 - x0) * (ai - 8.0) * temp; } } } /* objfun */ static void confun(Integer n, Integer ncnlin, Integer needc[], double x[], double conf[], double cjac[], Nag_Comm *comm) { #define CJAC(I,J) cjac[(I)*n + (J)] /* Function to evaluate the nonlinear constraints and its 1st derivatives. */ if (comm->first == Nag_TRUE) { /* First call to confun. Set all Jacobian elements to zero. * Note that this will only work when options.obj_deriv = TRUE * (the default). */ CJAC(0,0) = CJAC(0,1) = 0.0; } if (needc[0] > 0) { conf[0] = -0.09 - x[0]*x[1] + 0.49*x[1]; if (comm->flag == 2) { CJAC(0,0) = -x[1]; CJAC(0,1) = -x[0] + 0.49; } } } /* confun */ int main(void) { Integer exit_status_ex1=0; Integer exit_status_ex2=0; Vprintf("nag_opt_nlin_lsq (e04unc) Example Program Results\n"); exit_status_ex1 = ex1(); exit_status_ex2 = ex2(); return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1; } #undef FJAC #define A(I,J) a[(I)*tda + J] #define FJAC(I,J) fjac[(I)*tdfjac + J] static int ex1(void) { /* Local variables */ Integer exit_status=0, i, j, m, n, nbnd, nclin, ncnlin, tda, tdfjac; NagError fail; double *a=0, *bl=0, *bu=0, *f=0, *fjac=0, objf, *x=0, *y=0; INIT_FAIL(fail); Vprintf("\nExample 1: default options\n"); Vscanf(" %*[^\n]"); /* Skip heading in data file */ Vscanf(" %*[^\n]"); /* Read problem dimensions */ Vscanf(" %*[^\n]"); Vscanf("%ld%ld%*[^\n]", &m, &n); Vscanf(" %*[^\n]"); Vscanf("%ld%ld%*[^\n]", &nclin, &ncnlin); if (m>0 && n>0 && nclin>=0 && ncnlin>=0) { nbnd =n + nclin + ncnlin; if ( !( x = NAG_ALLOC(n, double)) || !( a = NAG_ALLOC(nclin*n, double)) || !( f = NAG_ALLOC(m, double)) || !( y = NAG_ALLOC(m, double)) || !( fjac = NAG_ALLOC(m*n, double)) || !( bl = NAG_ALLOC(nbnd, double)) || !( bu = NAG_ALLOC(nbnd, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdfjac = n; } else { Vprintf("Invalid m or n or nclin or ncnlin.\n"); exit_status = 1; return exit_status; } /* Read a, y, bl, bu and x from data file */ if (nclin > 0) { Vscanf(" %*[^\n]"); for (i = 0; i < nclin; ++i) for (j = 0; j < n; ++j) Vscanf("%lf",&A(i,j)); } /* Read the y vector of the objective */ Vscanf(" %*[^\n]"); for (i = 0; i < m; ++i) Vscanf("%lf",&y[i]); /* Read lower bounds */ Vscanf(" %*[^\n]"); for (i = 0; i < n + nclin + ncnlin; ++i) Vscanf("%lf",&bl[i]); /* Read upper bounds */ Vscanf(" %*[^\n]"); for (i = 0; i < n + nclin + ncnlin; ++i) Vscanf("%lf",&bu[i]); /* Read the initial point x */ Vscanf(" %*[^\n]"); for (i = 0; i < n; ++i) Vscanf("%lf",&x[i]); /* Solve the problem */ /* nag_opt_nlin_lsq (e04unc). * Solves nonlinear least-squares problems using the * sequential QP method */ nag_opt_nlin_lsq(m, n, nclin, ncnlin, a, tda, bl, bu, y, objfun, confun, x, &objf, f, fjac, tdfjac, E04_DEFAULT, NAGCOMM_NULL, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_nlin_lsq (e04unc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (x) NAG_FREE(x); if (a) NAG_FREE(a); if (f) NAG_FREE(f); if (y) NAG_FREE(y); if (fjac) NAG_FREE(fjac); if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); return exit_status; } /* ex1 */ static void user_print(const Nag_Search_State *st, Nag_Comm *comm) { static const char *states[] = {"IL", "IU", "FR", "LL", "UL", "EQ"}; Integer i, ixl, ixn; Nag_SimSt c_sim, f_sim; /* Derivative checking (only handle simple checks here ) */ if (comm->g_prt) { /* First ensure that a check has been performed */ if (st->gprint->g_chk.type != Nag_NoCheck) { Vprintf("\n\nDerivative Checks\n-----------------\n\n"); /* Objective check */ f_sim = st->gprint->f_sim; Vprintf("Simple check of objective gradients: "); if (f_sim.correct) Vprintf("OK\n"); else Vprintf("ERROR!\n"); Vprintf("Max error was %9.2e in subfunction %1ld\n\n", f_sim.max_error, f_sim.max_subfunction); /* Similarly for constraints */ c_sim = st->gprint->c_sim; Vprintf("Simple check of constraint gradients: "); if (c_sim.correct) Vprintf("OK\n"); else Vprintf("ERROR!\n"); Vprintf("Max error was %9.2e in constraint %1ld\n", c_sim.max_error, c_sim.max_constraint); } } /* Major iteration output */ if (comm->it_maj_prt) { if (st->first) { Vprintf("\nIterations\n----------\n"); Vprintf("\n Maj Mnr Step Merit function\n"); } Vprintf(" %4ld %4ld %8.1e %14.6e\n", st->iter, st->minor_iter, st->step, st->merit); } /* Solution output */ if (comm->sol_sqp_prt) { Vprintf("\nSolution\n--------\n"); /* Variable results */ Vprintf("\nVarbl State Value Lagr Mult\n"); for (i = 0; i < st->n; ++i) Vprintf("V%4ld %2s %14.6e %12.4e\n", i+1, states[st->state[i] + 2], st->x[i], st->lambda[i]); /* Linear constraint results */ Vprintf("\nL Con State Value Lagr Mult\n"); for (i = st->n; i < st->n + st->nclin; ++i) { ixl = i - st->n; Vprintf("V%4ld %2s %14.6e %12.4e\n", ixl+1, states[st->state[i] + 2], st->ax[ixl], st->lambda[i]); } /* Nonlinear constraint results */ Vprintf("\nN Con State Value Lagr Mult\n"); for (i = st->n + st->nclin; i < st->n + st->nclin + st->ncnlin; ++i) { ixn = i - st->n - st->nclin; Vprintf("V%4ld %2s %14.6e %12.4e\n", ixn+1, states[st->state[i] + 2], st->cx[ixn], st->lambda[i]); } } } /* user_print */ static int ex2(void) { /* Local variables */ Integer exit_status=0, i, j, m, n, nbnd, nclin, ncnlin, tda, tdfjac; NagError fail; Nag_E04_Opt options; double *a=0, *bl=0, *bu=0, *f=0, *fjac=0, objf, *x=0, *y=0; INIT_FAIL(fail); Vprintf("\nExample 2: user defined printing option\n"); Vscanf(" %*[^\n]"); /* Skip heading in data file */ /* Read problem dimensions */ Vscanf(" %*[^\n]"); Vscanf("%ld%ld%*[^\n]", &m, &n); Vscanf(" %*[^\n]"); Vscanf("%ld%ld%*[^\n]", &nclin, &ncnlin); if (m>0 && n>0 && nclin>=0 && ncnlin>=0) { nbnd = n + nclin + ncnlin; if ( !( x = NAG_ALLOC(n, double)) || !( a = NAG_ALLOC(nclin*n, double)) || !( f = NAG_ALLOC(m, double)) || !( y = NAG_ALLOC(m, double)) || !( fjac = NAG_ALLOC(m*n, double)) || !( bl = NAG_ALLOC(nbnd, double)) || !( bu = NAG_ALLOC(nbnd, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdfjac = n; } else { Vprintf("Invalid m or n nclin or ncnlin.\n"); exit_status = 1; return exit_status; } /* Read a, y, bl, bu and x from data file */ if (nclin > 0) { Vscanf(" %*[^\n]"); for (i = 0; i < nclin; ++i) for (j = 0; j < n; ++j) Vscanf("%lf",&A(i,j)); } /* Read the y vector of the objective */ Vscanf(" %*[^\n]"); for (i = 0; i < m; ++i) Vscanf("%lf",&y[i]); /* Read lower bounds */ Vscanf(" %*[^\n]"); for (i = 0; i < n + nclin + ncnlin; ++i) Vscanf("%lf",&bl[i]); /* Read upper bounds */ Vscanf(" %*[^\n]"); for (i = 0; i < n + nclin + ncnlin; ++i) Vscanf("%lf",&bu[i]); /* Read the initial point x */ Vscanf(" %*[^\n]"); for (i = 0; i < n; ++i) Vscanf("%lf",&x[i]); /* Set an option */ /* nag_opt_init (e04xxc). * Initialization function for option setting */ nag_opt_init(&options); options.print_fun = user_print; /* Solve the problem */ /* nag_opt_nlin_lsq (e04unc), see above. */ nag_opt_nlin_lsq(m, n, nclin, ncnlin, a, tda, bl, bu, y, objfun, confun, x, &objf, f, fjac, tdfjac, &options, NAGCOMM_NULL, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_nlin_lsq (e04unc).\n%s\n", fail.message); exit_status = 1; } /* nag_opt_free (e04xzc). * Memory freeing function for use with option setting */ nag_opt_free(&options, "all", &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (x) NAG_FREE(x); if (a) NAG_FREE(a); if (f) NAG_FREE(f); if (y) NAG_FREE(y); if (fjac) NAG_FREE(fjac); if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); return exit_status; } /* ex2 */