/* nag_opt_lsq_uncon_mod_func_comp (e04fcc) Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
*
* Mark 29.1, 2023.
*
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL lsqfun(Integer m, Integer n, const double x[],
double fvec[], Nag_Comm *comm);
static void NAG_CALL lsqgrd(Integer m, Integer n, double *fvec, double *fjac,
Integer ldfjac, double *g);
#ifdef __cplusplus
}
#endif
#define MMAX 15
#define TMAX 3
/* Define a user structure template to store data in lsqfun. */
struct user {
double y[MMAX];
double t[MMAX][TMAX];
};
int main(void) {
const char *optionsfile = "e04fcce.opt";
Integer exit_status = 0;
Nag_Boolean print;
Integer i, j, m, n, nt, tdfjac;
Nag_Comm comm;
Nag_E04_Opt options;
double *fjac = 0, fsumsq, *fvec = 0, *x = 0, *g = 0;
struct user s;
NagError fail;
INIT_FAIL(fail);
printf("nag_opt_lsq_uncon_mod_func_comp (e04fcc) Example Program Results\n");
fflush(stdout);
scanf(" %*[^\n]"); /* Skip heading in data file */
n = 3;
m = 15;
if (m >= 1 && n <= m) {
if (!(fjac = NAG_ALLOC(m * n, double)) || !(fvec = NAG_ALLOC(m, double)) ||
!(x = NAG_ALLOC(n, double)) || !(g = NAG_ALLOC(n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdfjac = n;
} else {
printf("Invalid m or n.\n");
exit_status = 1;
return exit_status;
}
/* Read data into structure.
* Observations t (j = 0, 1, 2) are held in s->t[i][j]
* (i = 0, 1, 2, . . ., 14)
*/
nt = 3;
for (i = 0; i < m; ++i) {
scanf("%lf", &s.y[i]);
for (j = 0; j < nt; ++j)
scanf("%lf", &s.t[i][j]);
}
/* Set up the starting point */
x[0] = 0.5;
x[1] = 1.0;
x[2] = 1.5;
/* nag_opt_init (e04xxc).
* Initialization function for option setting
*/
nag_opt_init(&options); /* Initialize options structure */
/* Set one option directly. */
/* nag_machine_precision (x02ajc).
* The machine precision
*/
options.optim_tol = 10.0 * sqrt(nag_machine_precision);
/* Read remaining option values from file */
print = Nag_FALSE;
/* nag_opt_read (e04xyc).
* Read options from a text file
*/
nag_opt_read("e04fcc", optionsfile, &options, print, "stdout", &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_read (e04xyc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Assign address of user defined structure to
* comm.p for communication to lsqfun().
*/
comm.p = (Pointer)&s;
/* nag_opt_lsq_uncon_mod_func_comp (e04fcc), see above. */
nag_opt_lsq_uncon_mod_func_comp(m, n, lsqfun, x, &fsumsq, fvec, fjac, tdfjac,
&options, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error/Warning from nag_opt_lsq_uncon_mod_func_comp (e04fcc).\n%s\n",
fail.message);
if (fail.code != NW_COND_MIN)
exit_status = 1;
}
if (fail.code == NE_NOERROR || fail.code == NW_COND_MIN) {
printf("On exit, the sum of squares is %12.4f\n", fsumsq);
printf("at the point");
for (i = 0; i < n; i++)
printf("%12.4f", x[i]);
printf("\n");
lsqgrd(m, n, fvec, fjac, tdfjac, g);
printf("The estimated gradient is");
for (i = 0; i < n; i++)
printf("%13.4e", g[i]);
printf("\n");
printf(" (machine dependent)\n");
printf("and the residuals are\n");
for (i = 0; i < m; i++)
printf("%9.1e\n", fvec[i]);
}
/* Free memory allocated to pointers s and v */
/* nag_opt_free (e04xzc).
* Memory freeing function for use with option setting
*/
nag_opt_free(&options, "all", &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
exit_status = 2;
goto END;
}
END:
NAG_FREE(fjac);
NAG_FREE(fvec);
NAG_FREE(x);
NAG_FREE(g);
return exit_status;
}
static void NAG_CALL lsqfun(Integer m, Integer n, const double x[],
double fvec[], Nag_Comm *comm) {
/* Function to evaluate the residuals.
*
* To avoid the use of a global varibale this example assigns the address
* of a user defined structure to comm.p in the main program (where the
* data was also read in).
* The address of this structure is recovered in each call to lsqfun()
* from comm->p and the structure used in the calculation of the residuals.
*/
Integer i;
struct user *s = (struct user *)comm->p;
for (i = 0; i < m; ++i)
fvec[i] =
x[0] + s->t[i][0] / (x[1] * s->t[i][1] + x[2] * s->t[i][2]) - s->y[i];
} /* lsqfun */
static void NAG_CALL lsqgrd(Integer m, Integer n, double *fvec, double *fjac,
Integer ldfjac, double *g) {
/* Function to evaluate gradient of the sum of squares */
NagError fail;
Integer i;
INIT_FAIL(fail);
nag_blast_dgemv(Nag_RowMajor, Nag_Trans, m, n, 1.0, fjac, ldfjac, fvec, 1,
0.0, g, 1, &fail);
for (i = 0; i < n; i++)
g[i] = 2.0 * g[i];
return;
}