/* nag_opt_handle_solve_bxnl (e04ggc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.6, 2022.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL lsqfun(Integer nvar, const double x[], Integer nres,
double rx[], Integer *inform, Nag_Comm *comm);
static void NAG_CALL lsqgrd(Integer nvar, const double x[], Integer nres,
Integer nnzrd, double rdx[], Integer *inform,
Nag_Comm *comm);
static void NAG_CALL lsqhes(Integer nvar, const double x[], Integer nres,
const double lambda[], double hx[], Integer *inform,
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
/* Model data */
const double t[24] = {0.00E+0, 5.00E-2, 1.00E-1, 1.50E-1, 2.00E-1, 2.50E-1,
3.00E-1, 3.50E-1, 4.00E-1, 4.50E-1, 5.00E-1, 5.50E-1,
6.00E-1, 6.50E-1, 7.00E-1, 7.50E-1, 8.00E-1, 8.50E-1,
9.00E-1, 9.50E-1, 1.00E+0, 1.05E+0, 1.10E+0, 1.15E+0};
const double y[24] = {2.5134, 2.0443, 1.6684, 1.3664, 1.1232, 0.9269,
0.7679, 0.6389, 0.5338, 0.4479, 0.3776, 0.3197,
0.2720, 0.2325, 0.1997, 0.1723, 0.1493, 0.1301,
0.1138, 0.1000, 0.0883, 0.0783, 0.0698, 0.0624};
int main(void) {
const double infbnd = 1.0e20;
Integer nvar, nres, isparse, nnzrd;
double x[6] = {1.2, 0.3, 5.6, 5.5, 6.5, 7.6};
double rinfo[100], stats[100];
double *rx, *blx, *bux, *z;
void *handle = 0;
Integer exit_status = 0;
/* Nag Types */
Nag_Comm comm;
NagError fail;
printf("nag_opt_handle_solve_bxnl (e04ggc) Example Program Results\n\n");
fflush(stdout);
nvar = 6;
nres = 24;
/* Allocate memory */
if (!(comm.user = NAG_ALLOC(2 * nres, double)) ||
!(blx = NAG_ALLOC(nvar, double)) || !(bux = NAG_ALLOC(nvar, double)) ||
!(rx = NAG_ALLOC(nres, double)) || !(z = NAG_ALLOC(nvar, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Fill the problem data structure
* copy t[]
*/
memcpy(&comm.user[0], t, 24 * sizeof(double));
/* copy y[] = */
memcpy(&comm.user[24], y, 24 * sizeof(double));
/* nag_opt_handle_init (e04rac).
* Initialize the handle
*/
nag_opt_handle_init(&handle, nvar, NAGERR_DEFAULT);
/* nag_opt_handle_set_nlnls (e04rmc)
* Define residuals structure, isparse=0 means the residual structure is
* dense => irowrd and icolrd arguments can be NULL
*/
isparse = 0;
nnzrd = 1;
nag_opt_handle_set_nlnls(handle, nres, isparse, nnzrd, NULL, NULL,
NAGERR_DEFAULT);
/* nag_opt_handle_opt_set (e04zmc)
* Set options
*/
/* Use second order information to build the TR model */
nag_opt_handle_opt_set(handle, "BXNL Use Second Derivatives = Yes",
NAGERR_DEFAULT);
/* Choose the Gauss-Newton method to solve the TR problem */
nag_opt_handle_opt_set(handle, "BXNL Model = Gauss-Newton", NAGERR_DEFAULT);
/* Add regularization */
nag_opt_handle_opt_set(handle, "BXNL Glob Method = Reg", NAGERR_DEFAULT);
/* Reduce printed output info */
nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);
/* Define bounds */
blx[0] = 0.0;
bux[0] = 1.0;
blx[1] = -1.0;
bux[1] = infbnd;
blx[2] = -1.0;
bux[2] = infbnd;
blx[3] = -1.0;
bux[3] = infbnd;
blx[4] = -1.0;
bux[4] = 1.0;
blx[5] = -1.0;
bux[5] = 10.0;
/* nag_opt_handle_set_simplebounds (e04rhc) */
nag_opt_handle_set_simplebounds(handle, nvar, blx, bux, NAGERR_DEFAULT);
/* nag_opt_handle_solve_bxnl (e04ggc)
* Call the solver
*/
INIT_FAIL(fail);
nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, lsqhes, NULLFN, NULLFN,
nvar, x, nres, rx, rinfo, stats, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_handle_solve_bxnl (e04ggc).\n%s\n",
fail.message);
exit_status = 1;
} else {
/* Retrieve the solution point */
nag_opt_handle_set_get_real(handle, "X", 1, &nvar, z, NAGERR_DEFAULT);
printf("\nSolver stored solution iterate in the handle\n");
printf("X: %8.2e %8.2e %8.2e %8.2e %8.2e %8.2e\n", z[0], z[1], z[2], z[3],
z[4], z[5]);
NAG_FREE(z);
}
END:
/* Free allocated memory */
if (handle)
/* nag_opt_handle_free (e04rzc).
* Free the problem handle and deallocate all the memory used
*/
nag_opt_handle_free(&handle, NAGERR_DEFAULT);
if (comm.user)
NAG_FREE(comm.user);
NAG_FREE(rx);
NAG_FREE(blx);
NAG_FREE(bux);
return exit_status;
}
static void NAG_CALL lsqfun(Integer nvar, const double x[], Integer nres,
double rx[], Integer *inform, Nag_Comm *comm) {
int i;
/* Interrupt the solver if the comm structure is not correctly initialized */
if (!comm || !(comm->user)) {
*inform = -1;
return;
}
/* Interrupt the solver if the data does not correspond to the problem */
if (nvar != 6 || nres != 24) {
*inform = -1;
return;
}
*inform = 0;
/* Fill the residuals values */
for (i = 0; i < nres; i++) {
rx[i] = comm->user[nres + i];
rx[i] = rx[i] - x[0] * exp(-x[1] * comm->user[i]) -
x[2] * exp(-x[3] * comm->user[i]) -
x[4] * exp(-x[5] * comm->user[i]);
}
}
static void NAG_CALL lsqgrd(Integer nvar, const double x[], Integer nres,
Integer nnzrd, double rdx[], Integer *inform,
Nag_Comm *comm) {
int i;
/* Interrupt the solver if the comm structure is not correctly initialized */
if (!comm || !(comm->user)) {
*inform = -1;
return;
}
/* Interrupt the solver if the data does not correspond to the problem */
if (nvar != 6 || nres != 24 || nnzrd != nres * nvar) {
*inform = -1;
return;
}
*inform = 0;
for (i = 0; i < nres; i++) {
rdx[i * nvar + 0] = -exp(-x[1] * comm->user[i]);
rdx[i * nvar + 1] = +comm->user[i] * x[0] * exp(-x[1] * comm->user[i]);
rdx[i * nvar + 2] = -exp(-x[3] * comm->user[i]);
rdx[i * nvar + 3] = +comm->user[i] * x[2] * exp(-x[3] * comm->user[i]);
rdx[i * nvar + 4] = -exp(-x[5] * comm->user[i]);
rdx[i * nvar + 5] = +comm->user[i] * x[4] * exp(-x[5] * comm->user[i]);
}
}
static void NAG_CALL lsqhes(Integer nvar, const double x[], Integer nres,
const double lambda[], double hx[], Integer *inform,
Nag_Comm *comm) {
double sum21 = 0.0, sum22 = 0.0, sum43 = 0.0, sum44 = 0.0, sum65 = 0.0,
sum66 = 0.0;
int i;
/* Interrupt the solver if the comm structure is not correctly initialized */
if (!comm || !(comm->user)) {
*inform = -1;
return;
}
/* Interrupt the solver if the data does not correspond to the problem */
if (nvar != 6 || nres != 24) {
*inform = -1;
return;
}
*inform = 0;
for (i = 0; i < nvar * nvar; i++)
hx[i] = 0.0;
for (i = 0; i < nres; i++) {
sum21 = sum21 + (lambda[i] * comm->user[i] * exp(-x[1] * comm->user[i]));
sum22 = sum22 + (-lambda[i] * pow(comm->user[i], 2) * x[0] *
exp(-x[1] * comm->user[i]));
sum43 = sum43 + (lambda[i] * comm->user[i] * exp(-x[3] * comm->user[i]));
sum44 = sum44 + (-lambda[i] * pow(comm->user[i], 2) * x[2] *
exp(-x[3] * comm->user[i]));
sum65 = sum65 + (lambda[i] * comm->user[i] * exp(-x[5] * comm->user[i]));
sum66 = sum66 + (-lambda[i] * pow(comm->user[i], 2) * x[4] *
exp(-x[5] * comm->user[i]));
}
hx[1 + 0 * nvar] = sum21;
hx[0 + 1 * nvar] = hx[1 + 0 * nvar];
hx[1 + 1 * nvar] = sum22;
hx[3 + 2 * nvar] = sum43;
hx[2 + 3 * nvar] = hx[3 + 2 * nvar];
hx[3 + 3 * nvar] = sum44;
hx[5 + 4 * nvar] = sum65;
hx[4 + 5 * nvar] = hx[5 + 4 * nvar];
hx[5 + 5 * nvar] = sum66;
}