/* nag_matop_real_nmf (f01sac) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */
#include <nag.h>

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

int main(void)
{
  /* Scalars */
  Integer exit_status = 0;
  Integer i, j, m, n, k, pda, pdw, pdh, seed, maxit;
  double errtol, norm, norma;
  /* Arrays */
  double *a = 0, *w = 0, *h = 0;
  /* Nag Types */
  Nag_OrderType order = Nag_ColMajor;
  NagError fail;

  INIT_FAIL(fail);

  /* Output preamble */
  printf("nag_matop_real_nmf (f01sac) ");
  printf("Example Program Results\n\n");
  fflush(stdout);

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

  /* Read in the matrix size and the required rank */
  scanf("%" NAG_IFMT "", &m);
  scanf("%" NAG_IFMT "", &n);
  scanf("%" NAG_IFMT "", &k);
  scanf("%*[^\n]");

  pda = m;
  if (!(a = NAG_ALLOC(pda * n, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  pdw = m;
  if (!(w = NAG_ALLOC(pdw * k, double)))
  {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }
  pdh = k;
  if (!(h = NAG_ALLOC(pda * n, double)))
  {
    printf("Allocation failure\n");
    exit_status = -3;
    goto END;
  }

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

  /* Choose the values of seed and errtol */
  seed = 23;
  errtol = 1.0e-6;

  /* Use the default value of maxit */
  maxit = -1;

  /* Find a non-negative matrix factorization using
   * nag_matop_real_nmf (f01sac) 
   * Non-negative matrix factorization of real non-negative matrix
   */
  nag_matop_real_nmf(m, n, k, a, pda, w, pdw, h, pdh, seed, errtol,
                     maxit, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_matop_real_nmf (f01sac)\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Print matrix W using nag_file_print_matrix_real_gen (x04cac)
   * Print real general matrix (easy-to-use)
   */
  nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag,
                         m, k, w, pdw, "W", NULL, &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;
  }

  printf("\n");

  /* Print matrix H using nag_file_print_matrix_real_gen (x04cac)
   * Print real general matrix (easy-to-use)
   */
  nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag,
                         k, n, h, pdh, "H", NULL, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen (x04cac)\n%s\n",
           fail.message);
    exit_status = 3;
    goto END;
  }
 
 /* Compute ||A||_F using
  * nag_blast_dge_norm (f16rac)
  */
  nag_blast_dge_norm(order, Nag_FrobeniusNorm, m, n, a, pda, &norma, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dge_norm (f16yac)\n%s\n", fail.message);
    exit_status = 4;
    goto END;
  }

 /* Compute residual A-WH using nag_blast_dgemm (f16yac)
  * Matrix-matrix product, two rectangular matrices
  */
  nag_blast_dgemm(order, Nag_NoTrans, Nag_NoTrans, m, n, k, -1.0, w, pdw, h, pdh,
            1.0, a, pda, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dgemm (f16yac)\n%s\n", fail.message);
    exit_status = 5;
    goto END;
  }

 /* Compute ||A-WH||_F using
  * nag_blast_dge_norm (f16rac)
  */
  nag_blast_dge_norm(order, Nag_FrobeniusNorm, m, n, a, pda, &norm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dge_norm (f16yac)\n%s\n", fail.message);
    exit_status = 6;
    goto END;
  }

  /* Print relative residual norm */
  printf("\nThe relative residual norm, ||A-WH||/||A||, is: %9.2e\n",
         norm/norma);

END:
  NAG_FREE(a);
  NAG_FREE(w);
  NAG_FREE(h);
  return exit_status;
}