/* nag_opt_lsq_check_deriv (e04yac) Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
*
* Mark 27.0, 2019.
*
*/
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL lsqfun(Integer m, Integer n, const double x[],
double fvec[], double fjac[], Integer tdfjac,
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
#define Y(I) comm.user[I]
#define T(I, J) comm.user[(I) *n + (J) + m]
#define YC(I) comm->user[(I)]
#define TC(I, J) comm->user[(I) *n + (J) + m]
#define FJAC(I, J) fjac[(I) *tdfjac + (J)]
int main(void)
{
Integer exit_status = 0, i, j, m, n, tdfjac;
NagError fail;
Nag_Comm comm;
double *fjac = 0, *fvec = 0, *work = 0, *x = 0;
INIT_FAIL(fail);
printf("nag_opt_lsq_check_deriv (e04yac) Example Program Results\n");
scanf(" %*[^\n]"); /* Skip heading in data file */
n = 3;
m = 15;
if (n >= 1 && m >= 1 && n <= m) {
if (!(fjac = NAG_ALLOC(m * n, double)) ||
!(fvec = NAG_ALLOC(m, double)) ||
!(x = NAG_ALLOC(n, double)) || !(work = NAG_ALLOC(m + m * n, double))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdfjac = n;
}
else {
printf("Invalid n or m.\n");
exit_status = 1;
return exit_status;
}
/* Allocate memory to communication array */
comm.user = work;
/* Observations t (j = 0, 1, 2) are held in T(i, j)
* (i = 0, 1, 2, . . ., 14) */
for (i = 0; i < m; ++i) {
scanf("%lf", &Y(i));
for (j = 0; j < n; ++j)
scanf("%lf", &T(i, j));
}
/* Set up an arbitrary point at which to check the 1st derivatives */
x[0] = 0.19;
x[1] = -1.34;
x[2] = 0.88;
printf("\nThe test point is ");
for (j = 0; j < n; ++j)
printf(" %12.3e", x[j]);
printf("\n");
/* nag_opt_lsq_check_deriv (e04yac).
* Least squares derivative checker for use with
* nag_opt_lsq_uncon_quasi_deriv_comp (e04gbc)
*/
nag_opt_lsq_check_deriv(m, n, lsqfun, x, fvec, fjac, tdfjac, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_lsq_check_deriv (e04yac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\nDerivatives are consistent with residual values.\n");
printf("\nAt the test point, lsqfun() gives\n\n");
printf(" Residuals 1st derivatives\n");
for (i = 0; i < m; ++i) {
printf(" %12.3e ", fvec[i]);
for (j = 0; j < n; ++j)
printf(" %12.3e", FJAC(i, j));
printf("\n");
}
END:
NAG_FREE(fjac);
NAG_FREE(fvec);
NAG_FREE(x);
NAG_FREE(work);
return exit_status;
}
static void NAG_CALL lsqfun(Integer m, Integer n, const double x[],
double fvec[], double fjac[], Integer tdfjac,
Nag_Comm *comm)
{
/* Function to evaluate the residuals and their 1st derivatives. */
Integer i;
double denom, dummy;
for (i = 0; i < m; ++i) {
denom = x[1] * TC(i, 1) + x[2] * TC(i, 2);
if (comm->flag != 1)
fvec[i] = x[0] + TC(i, 0) / denom - YC(i);
if (comm->flag != 0) {
FJAC(i, 0) = 1.0;
dummy = -1.0 / (denom * denom);
FJAC(i, 1) = TC(i, 0) * TC(i, 1) * dummy;
FJAC(i, 2) = TC(i, 0) * TC(i, 2) * dummy;
}
}
} /* lsqfun */