Example description
/* nag_matop_real_nmf_rcomm (f01sbc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */

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