/* nag_tsa_multi_autocorr_part (g13dbc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */

#include <nag.h>
#include <stdio.h>

int main(void) {
  /* Scalars */
  double v0;
  Integer exit_status, i1, i, j, j1, k, nk, nl, ns, nvp, pdc0, pddb;
  NagError fail;

  /* Arrays */
  double *c0 = 0, *c = 0, *d = 0, *db = 0, *p = 0, *v = 0, *w = 0, *wb = 0;

#define C(I, J, K) c[((K - 1) * ns + (J - 1)) * ns + I - 1]
#define D(I, J, K) d[((K - 1) * ns + (J - 1)) * ns + I - 1]
#define W(I, J, K) w[((K - 1) * ns + (J - 1)) * ns + I - 1]
#define WB(I, J, K) wb[((K - 1) * ns + (J - 1)) * ns + I - 1]

#ifdef NAG_COLUMN_MAJOR
#define C0(I, J) c0[(J - 1) * pdc0 + I - 1]
#define DB(I, J) db[(J - 1) * pddb + I - 1]
#else
#define C0(I, J) c0[(I - 1) * pdc0 + J - 1]
#define DB(I, J) db[(I - 1) * pddb + J - 1]
#endif

  INIT_FAIL(fail);

  exit_status = 0;

  printf("nag_tsa_multi_autocorr_part (g13dbc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Read series length, and numbers of lags */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &ns, &nl, &nk);

  if (ns > 0 && nl > 0 && nk > 0) {
    /* Allocate arrays */
    if (!(c0 = NAG_ALLOC(ns * ns, double)) ||
        !(c = NAG_ALLOC(ns * ns * nl, double)) ||
        !(d = NAG_ALLOC(ns * ns * nk, double)) ||
        !(db = NAG_ALLOC(ns * ns, double)) || !(p = NAG_ALLOC(nk, double)) ||
        !(v = NAG_ALLOC(nk, double)) ||
        !(w = NAG_ALLOC(ns * ns * nk, double)) ||
        !(wb = NAG_ALLOC(ns * ns * nk, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    pdc0 = ns;
    pddb = ns;

    /* Read autocovariances */
    for (i = 1; i <= ns; ++i) {
      for (j = 1; j <= ns; ++j)
        scanf("%lf", &C0(i, j));
    }
    scanf("%*[^\n] ");

    for (k = 1; k <= nl; ++k) {
      for (i = 1; i <= ns; ++i) {
        for (j = 1; j <= ns; ++j)
          scanf("%lf", &C(i, j, k));
      }
    }
    scanf("%*[^\n] ");

    /* Call routine to calculate multivariate partial
       autocorrelation function */

    /* nag_tsa_multi_autocorr_part (g13dbc).
     * Multivariate time series, multiple squared partial
     * autocorrelations
     */
    nag_tsa_multi_autocorr_part(c0, c, ns, nl, nk, p, &v0, v, d, db, w, wb,
                                &nvp, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_tsa_multi_autocorr_part (g13dbc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    if (fail.code == NE_NOERROR || fail.code == NE_POS_DEF) {
      printf("\n");
      printf("Number of valid parameters =%10" NAG_IFMT "\n", nvp);

      printf("\n");
      printf("Multivariate partial autocorrelations\n");

      for (i1 = 1; i1 <= nk; ++i1) {
        printf("%13.5f", p[i1 - 1]);
        if (i1 % 5 == 0 || i1 == nk)
          printf("\n");
      }

      printf("\n");
      printf("Zero lag predictor error variance determinant\n");
      printf("followed by error variance ratios\n");
      printf("%12.5f", v0);

      for (i1 = 1; i1 <= nk; ++i1) {
        printf("%13.5f", v[i1 - 1]);
        if (i1 % 5 == 0 || i1 == nk)
          printf("\n");
      }

      printf("\n");
      printf("Prediction error variances\n");
      printf("\n");

      for (k = 1; k <= nk; ++k) {
        printf("Lag =%5" NAG_IFMT "\n", k);
        for (i = 1; i <= ns; ++i) {
          for (j1 = 1; j1 <= ns; ++j1) {
            printf("%13.5f", D(i, j1, k));
            if (j1 % 5 == 0 || j1 == ns)
              printf("\n");
          }
        }
        printf("\n");
      }

      printf("Last backward prediction error variances\n");
      printf("\n");
      printf("Lag =%5" NAG_IFMT "\n", nvp);

      for (i = 1; i <= ns; ++i) {
        for (j1 = 1; j1 <= ns; ++j1) {
          printf("%13.5f", DB(i, j1));
          if (j1 % 5 == 0 || j1 == ns)
            printf("\n");
        }
      }

      printf("\n");
      printf("Prediction coefficients\n");
      printf("\n");

      for (k = 1; k <= nk; ++k) {
        printf("Lag =%5" NAG_IFMT "\n", k);
        for (i = 1; i <= ns; ++i) {
          for (j1 = 1; j1 <= ns; ++j1) {
            printf("%13.5f", W(i, j1, k));
            if (j1 % 5 == 0 || j1 == ns)
              printf("\n");
          }
        }
        printf("\n");
      }
      printf("Backward prediction coefficients\n");
      printf("\n");

      for (k = 1; k <= nk; ++k) {
        printf("Lag =%5" NAG_IFMT "\n", k);
        for (i = 1; i <= ns; ++i) {
          for (j1 = 1; j1 <= ns; ++j1) {
            printf("%13.5f", WB(i, j1, k));
            if (j1 % 5 == 0 || j1 == ns)
              printf("\n");
          }
        }
        printf("\n");
      }
    }
  }

END:
  NAG_FREE(c0);
  NAG_FREE(c);
  NAG_FREE(d);
  NAG_FREE(db);
  NAG_FREE(p);
  NAG_FREE(v);
  NAG_FREE(w);
  NAG_FREE(wb);

  return exit_status;
}