/* nag_opt_lsq_covariance (e04ycc) Example Program
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 * Mark 7 revised, 2001.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <math.h>
#include <nag_stdlib.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C" {
#endif
static void lsqfun(Integer m, Integer n, double x[], double fvec[],
                   Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

#define MMAX 15
#define NMAX 3
#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)
{
  double fjac[MMAX][NMAX], fvec[MMAX], x[NMAX], cj[NMAX];
  Integer i, j, m, n, nt, tdj, job;
  double fsumsq;
  Nag_E04_Opt options;
  Nag_Comm comm;
  static NagError fail, fail2;
  struct user s;

  Vprintf("e04ycc Example Program Results.\n");
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  n = 3;
  m = 15;
  tdj = NMAX;
  nt = 3;

  /* 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)
    {
      Vscanf("%lf", &s.y[i]);
      for (j = 0; j < nt; ++j) Vscanf("%lf", &s.t[i][j]);
    }

  /* Set up the starting point */
  x[0] = 0.5;
  x[1] = 1.0;
  x[2] = 1.5;

  e04xxc(&options); /* Initialise options structure */

  /* Assign address of user defined structure to
   * comm.p for communication to lsqfun().
   */
  comm.p = (Pointer)&s;

  fail.print = TRUE;
  e04fcc(m, n, lsqfun, x, &fsumsq, fvec, &fjac[0][0], tdj,
         &options, &comm, &fail);

  if (fail.code == NE_NOERROR || fail.code == NW_COND_MIN)
    {
      job = 0;
      e04ycc(job, m, n, fsumsq, cj, &options, &fail);

      if (fail.code == NE_NOERROR)
        {
          Vprintf("\nEstimates of the variances of the sample regression");
          Vprintf(" coefficients are:\n");
          for (i = 0; i < n; ++i)
            Vprintf(" %15.5e", cj[i]);
          Vprintf("\n");
        }
    }

  /* Free memory allocated to pointers s and v */
  fail2.print = TRUE;
  e04xzc(&options, "all", &fail2);
  if (fail.code != NE_NOERROR || fail2.code != NE_NOERROR) return EXIT_FAILURE;
  return EXIT_SUCCESS;
}

static void lsqfun(Integer m, Integer n, double x[], double fvec[],
                   Nag_Comm *comm)
{
  /* Function to evaluate the residuals.
   *
   * The address of the user defined 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 */