NAG Library Manual, Mark 29.3
Interfaces:  FL   CL   CPP   AD 

NAG CL Interface Introduction
Example description
/* nag_matop_real_nmf_rcomm (f01sbc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */

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

#define A_DENSE(I, J) a_dense[(J - 1) * pda + I - 1]

int main(void) {
  Integer exit_status = 0, i, j, irevcm = 0, m, n, k, icount = 0, maxit,
          exit_loop = 0, pdw, pdh, pdht, pda, q, seed, nnz;
  double *comm = 0, *a_dense = 0, *a = 0, *w = 0, *h = 0, *ht = 0;
  Integer icomm[9], *icolzp = 0, *irowix = 0, element;
  double errtol, norm, norma;
  /* Nag Types */
  Nag_OrderType order = Nag_ColMajor;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_matop_real_nmf_rcomm (f01sbc) Example Program Results\n");

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

  /* Read values of m, n and k */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &m, &n, &k);

  /* Pretend A is a square q x q matrix to allow use of
   * nag_sparse_direct_real_gen_matmul for matrix multiplications
   */
  q = MAX(m, n);

  /* Allocate memory */
  pdw = q;
  pdh = k;
  pdht = q;
  pda = m;
  if (!(w = NAG_ALLOC(pdw * k, double)) || !(h = NAG_ALLOC(pdh * n, double)) ||
      !(ht = NAG_ALLOC(pdht * k, double)) ||
      !(a_dense = NAG_ALLOC(pda * n, double)) ||
      !(icolzp = NAG_ALLOC(q + 1, Integer)) ||
      !(comm = NAG_ALLOC((2 * m + n) * k + 3, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read the matrix A */
  for (i = 0; i < n + 1; ++i)
    scanf("%" NAG_IFMT "%*[^\n] ", &icolzp[i]);
  nnz = icolzp[n] - 1;

  if (!(a = NAG_ALLOC(nnz, double)) || !(irowix = NAG_ALLOC(nnz, Integer))) {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }
  for (i = 0; i < nnz; ++i)
    scanf("%lf%" NAG_IFMT "%*[^\n] ", &a[i], &irowix[i]);

  /* If m > n (q==m), nag_sparse_direct_real_gen_matmul will need
   * a total of q columns, otherwise no adjustment is necessary
   */
  for (i = n + 1; i < q + 1; ++i)
    icolzp[i] = icolzp[n];

  /* Initialize w and ht to zero */
  for (i = 0; i < pdw * k; ++i)
    w[i] = 0.0;

  for (i = 0; i < pdht * k; ++i)
    ht[i] = 0.0;

  /* Choose the values of seed and errtol */
  seed = 234;
  errtol = 1.0e-3;

  /* We will terminate after 200 iterations */
  maxit = 200;

  /* nag_matop_real_nmf_rcomm (f01sbc).
   * Non-negative matrix factorisation
   * (reverse communication)
   */
  do {
    nag_matop_real_nmf_rcomm(&irevcm, m, n, k, w, pdw, h, pdh, ht, pdht, seed,
                             errtol, comm, icomm, &fail);

    switch (irevcm) {
    case 1:
      if (icount == maxit) {
        printf("Exiting after the maximum number of iterations\n\n");
        exit_loop = 1;
      } else {
        icount++;
      }
      /* w and h are available for printing */
      break;
    case 2:
      /* compute a^t * w and store in ht */
      nag_sparse_direct_real_gen_matmul(order, Nag_Trans, q, k, 1.0, icolzp,
                                        irowix, a, w, pdw, 0.0, ht, pdht,
                                        &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_sparse_direct_real_gen_matmul (f11mkc) \
                   \n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
      break;
    case 3:
      /* compute a * ht and store in w */
      nag_sparse_direct_real_gen_matmul(order, Nag_NoTrans, q, k, 1.0, icolzp,
                                        irowix, a, ht, pdht, 0.0, w, pdw,
                                        &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_sparse_direct_real_gen_matmul (f11mkc) \
                   \n%s\n",
               fail.message);
        exit_status = 2;
        goto END;
      }
      break;
    }
  } while (irevcm != 0 && exit_loop == 0);

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_matop_real_nmf_rcomm (f01sbc).\n%s\n", fail.message);
    exit_status = 3;
    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 = 4;
    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 = 5;
    goto END;
  }

  /* In order to compute the residual we convert a to dense format */
  for (i = 1; i <= n; i++) {
    for (j = 1; j <= m; j++) {
      A_DENSE(j, i) = 0.0;
    }
  }

  element = 0;
  for (i = 1; i <= n; i++) {
    for (j = 1; j <= icolzp[i] - icolzp[i - 1]; j++) {
      A_DENSE(irowix[element], i) = a[element];
      element++;
    }
  }

  /* Compute ||A||_F using
   * nag_blast_dge_norm (f16rac)
   */
  nag_blast_dge_norm(order, Nag_FrobeniusNorm, m, n, a_dense, pda, &norma,
                     &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dge_norm (f16yac)\n%s\n", fail.message);
    exit_status = 6;
    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_dense, pda, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dgemm (f16yac)\n%s\n", fail.message);
    exit_status = 7;
    goto END;
  }

  /* Compute ||A-WH||_F using
   * nag_blast_dge_norm (f16rac)
   */
  nag_blast_dge_norm(order, Nag_FrobeniusNorm, m, n, a_dense, pda, &norm,
                     &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dge_norm (f16yac)\n%s\n", fail.message);
    exit_status = 8;
    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(a_dense);
  NAG_FREE(icolzp);
  NAG_FREE(irowix);
  NAG_FREE(comm);
  NAG_FREE(w);
  NAG_FREE(h);
  NAG_FREE(ht);

  return exit_status;
}