NAG Library Manual, Mark 28.4
/* nag_opt_lsq_check_deriv (e04yac) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.4, 2022.
*
*/

#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 */