/* nag_correg_pls_fit (g02lcc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */
/* Pre-processor includes */
#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>

int main(void) {
  /*Integer scalar and array declarations */
  Integer exit_status = 0;
  Integer i, ip, ip1, j, maxfac, my, nfact, vipopt;
  Integer pdb, pdc, pdob, pdp, pdvip, pdw, pdycv;
  /*Double scalar and array declarations */
  double rcond;
  double *b = 0, *c = 0, *ob = 0, *p = 0, *vip = 0, *w = 0;
  double *xbar = 0, *xstd = 0, *ybar = 0, *ycv = 0, *ystd = 0;
  /*Character scalar and array declarations */
  char siscale[40], sorig[40];
  /*NAG Types */
  Nag_OrderType order;
  Nag_ScalePredictor iscale;
  Nag_EstimatesOption orig;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_correg_pls_fit (g02lcc) Example Program Results\n");
  /* Skip header in data file */
  scanf("%*[^\n] ");
  /* Read data values */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%39s "
        "%39s %" NAG_IFMT "%*[^\n] ",
        &ip, &my, &maxfac, &nfact, sorig, siscale, &vipopt);
  orig = (Nag_EstimatesOption)nag_enum_name_to_value(sorig);
  iscale = (Nag_ScalePredictor)nag_enum_name_to_value(siscale);

#ifdef NAG_COLUMN_MAJOR
  pdb = ip;
  pdc = my;
#define C(I, J) c[(J - 1) * pdc + I - 1]
  pdob = ip + 1;
  pdp = ip;
#define P(I, J) p[(J - 1) * pdp + I - 1]
  pdvip = ip;
  pdw = ip;
#define W(I, J) w[(J - 1) * pdw + I - 1]
  pdycv = maxfac;
#define YCV(I, J) ycv[(J - 1) * pdycv + I - 1]
  order = Nag_ColMajor;
#else
  pdb = my;
  pdc = maxfac;
#define C(I, J) c[(I - 1) * pdc + J - 1]
  pdob = my;
  pdp = maxfac;
#define P(I, J) p[(I - 1) * pdp + J - 1]
  pdvip = vipopt;
  pdw = maxfac;
#define W(I, J) w[(I - 1) * pdw + J - 1]
  pdycv = my;
#define YCV(I, J) ycv[(I - 1) * pdycv + J - 1]
  order = Nag_RowMajor;
#endif
  if (!(b = NAG_ALLOC(pdb * (order == Nag_RowMajor ? ip : my), double)) ||
      !(c = NAG_ALLOC(pdc * (order == Nag_RowMajor ? my : maxfac), double)) ||
      !(ob = NAG_ALLOC(pdob * (order == Nag_RowMajor ? (ip + 1) : my),
                       double)) ||
      !(p = NAG_ALLOC(pdp * (order == Nag_RowMajor ? ip : maxfac), double)) ||
      !(vip =
            NAG_ALLOC(pdvip * (order == Nag_RowMajor ? ip : vipopt), double)) ||
      !(w = NAG_ALLOC(pdw * (order == Nag_RowMajor ? ip : maxfac), double)) ||
      !(xbar = NAG_ALLOC(ip, double)) || !(xstd = NAG_ALLOC(ip, double)) ||
      !(ybar = NAG_ALLOC(my, double)) ||
      !(ycv =
            NAG_ALLOC(pdycv * (order == Nag_RowMajor ? maxfac : my), double)) ||
      !(ystd = NAG_ALLOC(my, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Read P */
  for (i = 1; i <= ip; i++) {
    for (j = 1; j <= maxfac; j++)
      scanf("%lf ", &P(i, j));
  }
  scanf("%*[^\n] ");
  /* Read C */
  for (i = 1; i <= my; i++) {
    for (j = 1; j <= maxfac; j++)
      scanf("%lf ", &C(i, j));
  }
  scanf("%*[^\n] ");
  /* Read W */
  for (i = 1; i <= ip; i++) {
    for (j = 1; j <= maxfac; j++)
      scanf("%lf ", &W(i, j));
  }
  scanf("%*[^\n] ");
  /* Read YCV */
  for (i = 1; i <= maxfac; i++) {
    for (j = 1; j <= my; j++)
      scanf("%lf ", &YCV(i, j));
  }
  scanf("%*[^\n] ");
  /* Read means */
  if (orig == Nag_EstimatesOrig) {
    for (j = 0; j < ip; j++)
      scanf("%lf ", &xbar[j]);
    scanf("%*[^\n] ");
    for (j = 0; j < my; j++)
      scanf("%lf ", &ybar[j]);
    scanf("%*[^\n] ");
    if (iscale != Nag_PredNoScale) {
      for (j = 0; j < ip; j++)
        scanf("%lf ", &xstd[j]);
      scanf("%*[^\n] ");
      for (j = 0; j < my; j++)
        scanf("%lf ", &ystd[j]);
      scanf("%*[^\n] ");
    }
  }
  /* Calculate predictions */
  rcond = -1.00e0;
  ip1 = ip + 1;
  /*
   * nag_correg_pls_fit (g02lcc)
   * Partial least squares
   */
  nag_correg_pls_fit(order, ip, my, maxfac, nfact, p, pdp, c, pdc, w, pdw,
                     rcond, b, pdb, orig, xbar, ybar, iscale, xstd, ystd, ob,
                     pdob, vipopt, ycv, pdycv, vip, pdvip, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_correg_pls_fit (g02lcc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  /*
   * nag_file_print_matrix_real_gen (x04cac)
   * Print real general matrix (easy-to-use)
   */
  fflush(stdout);
  nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, ip,
                                 my, b, pdb, "B ", 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  if (orig == Nag_EstimatesOrig) {
    /*
     * nag_file_print_matrix_real_gen (x04cac)
     * Print real general matrix (easy-to-use)
     */
    fflush(stdout);
    nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag,
                                   ip1, my, ob, pdob, "OB", 0, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }
  if (vipopt != 0) {
    /*
     * nag_file_print_matrix_real_gen (x04cac)
     * Print real general matrix (easy-to-use)
     */
    fflush(stdout);
    nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag,
                                   ip, vipopt, vip, pdvip, "VIP", 0, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }

END:
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(ob);
  NAG_FREE(p);
  NAG_FREE(vip);
  NAG_FREE(w);
  NAG_FREE(xbar);
  NAG_FREE(xstd);
  NAG_FREE(ybar);
  NAG_FREE(ycv);
  NAG_FREE(ystd);

  return exit_status;
}