/* nag_rnla_randproj_dct_real (f10dac) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.7, 2022.
 */

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

int main(void) {

  /* Scalars */
  Integer exit_status = 0, i, j, m, n, pda, lda;
  /* Arrays */
  double *a = 0, *a_copy = 0;

  /* Nag Types */
  NagError fail;
  Nag_OrderType order;

  /* Declarations for random projection */
  Integer k, *d = 0, *col_ind = 0;
  double *pa = 0;
  char nag_enum_arg1[40];
  char nag_enum_arg2[40];
  Nag_TransType transa;

  /* Declarations for RNG */
  Integer lstate;
  Integer *state = 0;
  double *x = 0;
  Nag_BaseRNG genid = Nag_Basic;
  Integer subid = 0;
  Integer seed[] = {1762543};
  Integer lseed = 1;

  /* Declarations for lapack routines    */
  Integer tau_len, pdu, pdvt;
  double *tau = 0, *q = 0, *dgemm_inout = 0, *res_mat = 0, *work = 0, *s = 0,
         *u = 0, *vt = 0;
  double alpha, beta, sum;

  /* Declarations for error estimation   */
  Integer r;
  double xmu, var, err_bound, max_res_norm;
  double *omega = 0, *y = 0, *res = 0, *res_norm = 0;

#define A(I, J) a[(J - 1) * pda + I - 1]
  order = Nag_ColMajor;

  INIT_FAIL(fail);

  printf("nag_rnla_randproj_dct_real (f10dac) Example Program Results\n\n");
  /* Skip heading in data file */
  scanf("%*[^\n]");
  /* Read the problem dimensions, m, n, k */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &m, &n);

  if (m < 0 || n < 0) {
    printf("Invalid m or n\n");
    exit_status = 1;
    goto END;
  }

  /* Read in random projection routine parameters - mode, transpose, k */
  scanf("%39s %39s %" NAG_IFMT "%*[^\n] ", nag_enum_arg1, nag_enum_arg2, &k);

  /* nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  transa = (Nag_TransType)nag_enum_name_to_value(nag_enum_arg2);

  /* Allocate memory */
  if (!(a = NAG_ALLOC(m * n, double)) || !(a_copy = NAG_ALLOC(m * n, double)) ||
      !(pa = NAG_ALLOC(m * k, double)) || !(d = NAG_ALLOC(n, Integer)) ||
      !(col_ind = NAG_ALLOC(k, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

#ifdef NAG_COLUMN_MAJOR
  pda = m;
  pdu = m;
  pdvt = n;
#else
  pda = n;
  pdu = n;
  pdvt = m;
#endif

  /* Read the m by n matrix A from data file */
  for (i = 1; i <= m; ++i)
    for (j = 1; j <= n; ++j)
      scanf("%lf", &A(i, j));
  scanf("%*[^\n]");

  /* Copy a into a_copy */
  for (i = 0; i < m * n; i++)
    a_copy[i] = a[i];

  nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, m,
                                 n, a, pda, "Matrix A", 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
           fail.message);
    exit_status = 2;
    goto END;
  }

  /* Random number generation */

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

  /* Get the length of the state array */
  lstate = -1;
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 3;
    goto END;
  }

  /* Allocate arrays */
  if (!(x = NAG_ALLOC(n, double)) || !(state = NAG_ALLOC(lstate, Integer))) {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }

  /* Initialize the generator to a repeatable sequence */
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 4;
    goto END;
  }
  lda = m;
  /* Perform random projection */
  /* Copy a into a_copy */
  for (i = 0; i < m * n; i++)
    a_copy[i] = a[i];
  nag_rnla_randproj_dct_real(transa, m, n, a_copy, lda, k, state, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rnla_randproj_dct_real (f10dac).\n%s\n",
           fail.message);
    exit_status = 5;
    goto END;
  }
  for (i = 0; i < m * k; i++)
    pa[i] = a[i];

  /* Perform QR decomposition on random projection */
  tau_len = MIN(m, k);

  /* Allocate memory */
  if (!(tau = NAG_ALLOC(tau_len, double))) {
    printf("Allocation failure\n");
    exit_status = -3;
    goto END;
  }

  nag_lapackeig_dgeqrf(order, m, k, pa, m, tau, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapackeig_dgeqrf (f08aec).\n%s\n", fail.message);
    exit_status = 6;
    goto END;
  }
  q = pa;

  nag_lapackeig_dorgqr(order, m, k, k, pa, m, tau, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapackeig_dorgqr (f08afc).\n%s\n", fail.message);
    exit_status = 7;
    goto END;
  }

  /* Calculate true error as largest singular value of (I-QQ')A */

  /* Calculate (I - Q Q') */
  if (!(dgemm_inout = NAG_ALLOC(m * m, double)) ||
      !(res_mat = NAG_ALLOC(m * n, double)) ||
      !(s = NAG_ALLOC(MIN(m, n), double)) || !(u = NAG_ALLOC(m * m, double)) ||
      !(vt = NAG_ALLOC(m * n, double)) ||
      !(work = NAG_ALLOC(MIN(m, n), double))) {
    printf("Allocation failure\n");
    exit_status = -4;
    goto END;
  }

  /* Define identity matrix */
  for (i = 0; i < m; ++i) {
    for (j = 0; j < n; ++j) {
      dgemm_inout[i * m + j] = 0.0;
    }
    dgemm_inout[i * m + i] = 1.0;
  }

  alpha = -1.0;
  beta = 1.0;
  nag_blast_dgemm(order, Nag_NoTrans, Nag_Trans, m, m, k, alpha, q, m, q, m,
                  beta, dgemm_inout, m, &fail);

  /* Calculate (I - Q Q') A */
  alpha = 1.0;
  beta = 0.0;
  nag_blast_dgemm(order, Nag_NoTrans, Nag_NoTrans, m, n, m, alpha, dgemm_inout,
                  m, a, m, beta, res_mat, m, &fail);

  /* Compute the singular values of res_mat */
  nag_lapackeig_dgesvd(order, Nag_NotU, Nag_NotVT, m, n, res_mat, pda, s, u,
                       pdu, vt, pdvt, work, &fail);

  /* Calculate error estimate */
  r = 10;
  if (!(omega = NAG_ALLOC(n, double)) || !(y = NAG_ALLOC(m, double)) ||
      !(res = NAG_ALLOC(m, double)) || !(res_norm = NAG_ALLOC(r, double))) {
    printf("Allocation failure\n");
    exit_status = -5;
    goto END;
  }
  max_res_norm = 0;
  for (j = 0; j < r; ++j) {
    /* Generate Gaussian random vector, omega */
    xmu = 0.0;
    var = 1.0;
    nag_rand_dist_normal(n, xmu, var, state, omega, &fail);
    /* calculate y = A omega */
    nag_blast_dgemv(order, Nag_NoTrans, m, n, alpha, a, m, omega, 1, beta, y, 1,
                    &fail);
    /* Calculate (I - Q Q') y */
    nag_blast_dgemv(order, Nag_NoTrans, m, m, alpha, dgemm_inout, m, y, 1, beta,
                    res, 1, &fail);
    sum = 0;
    for (i = 0; i < m; ++i)
      sum = sum + pow(res[i], 2);

    res_norm[j] = sqrt(sum);
    if (res_norm[j] > max_res_norm) {
      max_res_norm = res_norm[j];
    }
  }

  err_bound = sqrt(200.0 / nag_math_pi) * max_res_norm;

  printf("\n True error \n %10.1e\n", s[0]);
  printf("\n Probabilistic error bound \n %10.1e\n", err_bound);
  printf("\n Maximum probability that true error exceeds error bound \n");
  printf(" %10.1e\n", pow(10.0, -r));
END:

  NAG_FREE(a);
  NAG_FREE(a_copy);
  NAG_FREE(pa);
  NAG_FREE(d);
  NAG_FREE(col_ind);
  NAG_FREE(x);
  NAG_FREE(state);
  NAG_FREE(tau);
  NAG_FREE(dgemm_inout);
  NAG_FREE(res_mat);
  NAG_FREE(s);
  NAG_FREE(u);
  NAG_FREE(vt);
  NAG_FREE(work);
  NAG_FREE(omega);
  NAG_FREE(y);
  NAG_FREE(res);
  NAG_FREE(res_norm);

  return exit_status;
}