/* nag_stat_pdf_multi_normal_vector (g01lbc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.0, 2023.
 */

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

#define X(I, J) x[(J)*pdx + I]
#define SIG(I, J) sig[(J)*pdsig + I]

int main(void) {
  /* Integer scalar and array declarations */
  Integer i, j, k, n, pdx, pdsig, rank;
  Integer exit_status = 0;

  /* NAG structures and types */
  NagError fail;
  Nag_MatrixType iuld;
  Nag_Boolean ilog;

  /* Double scalar and array declarations */
  double *x = 0, *xmu = 0, *sig = 0, *pdf = 0;

  /* Character scalar and array declarations */
  char ciuld[40], cilog[40];

  /* Initialize the error structure */
  INIT_FAIL(fail);

  printf(
      "nag_stat_pdf_multi_normal_vector (g01lbc) Example Program Results\n\n");
  fflush(stdout);

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

  /* Read in the problem size */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%39s%39s%*[^\n] ", &k, &n, ciuld, cilog);
  ilog = (Nag_Boolean)nag_enum_name_to_value(cilog);
  iuld = (Nag_MatrixType)nag_enum_name_to_value(ciuld);

  pdx = n;
  pdsig = (iuld == Nag_DiagonalMatrix) ? 1 : n;
  if (!(x = NAG_ALLOC(pdx * k, double)) || !(xmu = NAG_ALLOC(n, double)) ||
      !(sig = NAG_ALLOC(pdsig * n, double)) || !(pdf = NAG_ALLOC(k, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in and echo the vector of means */
  for (i = 0; i < n; i++)
    scanf("%lf", &xmu[i]);
  scanf("%*[^\n] ");

  printf("Vector of Means:\n");
  for (i = 0; i < n; i++)
    printf("%8.4f ", xmu[i]);
  printf("\n\n");
  fflush(stdout);

  /* Read in and echo the covariance matrix */

  if (iuld == Nag_DiagonalMatrix) {

    for (i = 0; i < n; i++)
      scanf("%lf", &SIG(1, i));
    scanf("%*[^\n] ");

    printf("Diagonal Elements of Covariance Matrix:\n");
    for (i = 0; i < n; i++)
      printf("%8.4f ", SIG(1, i));
    printf("\n\n");

  } else if (iuld == Nag_LowerMatrix || iuld == Nag_LowerFactored) {

    for (i = 0; i < n; i++) {
      for (j = 0; j <= i; j++)
        scanf("%lf", &SIG(i, j));
      scanf("%*[^\n] ");
    }
    /* Print details of Covariance matrix using
     * nag_file_print_matrix_real_gen (x04cac)
     */
    if (iuld == Nag_LowerMatrix)
      nag_file_print_matrix_real_gen(Nag_ColMajor, Nag_LowerMatrix,
                                     Nag_NonUnitDiag, n, n, sig, pdsig,
                                     "Covariance Matrix:", NULL, &fail);
    else
      nag_file_print_matrix_real_gen(Nag_ColMajor, Nag_LowerMatrix,
                                     Nag_NonUnitDiag, n, n, sig, pdsig,
                                     "Lower Triangular Cholesky Factor "
                                     "of Covariance Matrix:",
                                     NULL, &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;
    }

  } else {
    /* Upper triangle of matrix or factor is stored */

    for (i = 0; i < n; i++) {
      for (j = i; j < n; j++)
        scanf("%lf", &SIG(i, j));
      scanf("%*[^\n] ");
    }
    /* Print details of Covariance matrix using
     * nag_file_print_matrix_real_gen (x04cac)
     */
    if (iuld == Nag_UpperMatrix)
      nag_file_print_matrix_real_gen(Nag_ColMajor, Nag_UpperMatrix,
                                     Nag_NonUnitDiag, n, n, sig, pdsig,
                                     "Covariance Matrix:", NULL, &fail);
    else
      nag_file_print_matrix_real_gen(Nag_ColMajor, Nag_UpperMatrix,
                                     Nag_NonUnitDiag, n, n, sig, pdsig,
                                     "Upper Triangular Cholesky Factor "
                                     "of Covariance Matrix:",
                                     NULL, &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;
    }
  }

  /* Read in the points at which to evaluate the PDF */
  for (j = 0; j < k; j++)
    for (i = 0; i < n; i++)
      scanf("%lf", &X(i, j));

  /* nag_stat_pdf_multi_normal_vector (g01lbc): evaluate the PDF */
  nag_stat_pdf_multi_normal_vector(ilog, k, n, x, pdx, xmu, iuld, sig, pdsig,
                                   pdf, &rank, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_stat_pdf_multi_normal_vector (g01lbc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Display results */
  printf("\nRank of the covariance matrix: %" NAG_IFMT "\n\n", rank);
  if (ilog)
    printf("     log(PDF)                  X\n");
  else
    printf("       PDF                     X\n");
  printf("  ------------------------------------------------\n");
  for (j = 0; j < k; j++) {
    printf(" %13.4e", pdf[j]);
    for (i = 0; i < n; i++)
      printf(" %8.4f", X(i, j));
    printf("\n");
  }

END:
  NAG_FREE(x);
  NAG_FREE(xmu);
  NAG_FREE(sig);
  NAG_FREE(pdf);

  return exit_status;
}