/* nag_rnla_randproj_dct_real (f10dac) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */

#include <math.h>
#include <stdio.h>
#include <nag.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;
}