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

NAG CL Interface Introduction
Example description
/* nag_sparse_real_gen_sort (f11zac) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.0, 2024.
 */

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

int main(void) {
  double *b = 0, *c = 0, *a = 0, anorm = 0.0;
  Integer *bcol = 0, *ccol = 0, *icol = 0, *icol2 = 0, *trcol = 0;
  Integer *brow = 0, *crow = 0, *irow = 0, *istr = 0, *trrow = 0;
  Integer exit_status = 0, i, j, n, bnnz, cnnz, nnz;
  Nag_SparseNsym_Zeros zero;
  Nag_SparseNsym_Dups dup;
  Nag_NormType norm;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_sparse_real_gen_sort (f11zac) Example Program Results\n\n");

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

  /* Read order of matrices and number of nonzero entries of B */
  scanf("%" NAG_IFMT "%*[^\n]", &n);
  scanf("%" NAG_IFMT "%*[^\n]", &bnnz);

  /* Allocate memory */
  b = NAG_ALLOC(bnnz, double);
  brow = NAG_ALLOC(bnnz, Integer);
  bcol = NAG_ALLOC(bnnz, Integer);

  if (!b || !brow || !bcol) {
    printf("Allocation failure for b\n");
    exit_status = -1;
    goto END;
  }

  /* Read the nonzero elements of B */
  for (i = 0; i < bnnz; ++i)
    scanf("%lf%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &b[i], &brow[i], &bcol[i]);

    /* Read number of nonzero entries of C */
  scanf("%" NAG_IFMT "%*[^\n]", &cnnz);

  /* Allocate memory */
  c = NAG_ALLOC(cnnz, double);
  crow = NAG_ALLOC(cnnz, Integer);
  ccol = NAG_ALLOC(cnnz, Integer);

  if (!c || !crow || !ccol) {
    printf("Allocation failure for c\n");
    exit_status = -1;
    goto END;
  }

  /* Read the nonzero elements of C */
  for (i = 0; i < cnnz; ++i)
    scanf("%lf%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &c[i], &crow[i], &ccol[i]);

  /* Output the original non-zero elements */
  printf("Elements of B\n");
  printf("nnz =  %4" NAG_IFMT "\n", bnnz);
  printf("%8s%16s%8s%8s\n", "", "b", "brow", "bcol");
  for (i = 0; i < bnnz; ++i)
    printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "%8" NAG_IFMT "\n", i, b[i],
           brow[i], bcol[i]);

  printf("Elements of C\n");
  printf("nnz =  %4" NAG_IFMT "\n", cnnz);
  printf("%8s%16s%8s%8s\n", "", "c", "crow", "ccol");
  for (i = 0; i < cnnz; ++i)
    printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "%8" NAG_IFMT "\n", i, c[i],
           crow[i], ccol[i]);

  /* Set initial size for resulting array a, and allocate ememory */
  nnz = bnnz + cnnz;

  a = NAG_ALLOC(nnz, double);
  irow = NAG_ALLOC(nnz, Integer);
  icol = NAG_ALLOC(nnz, Integer);
  istr = NAG_ALLOC(n + 1, Integer);

  if (!a || !irow || !icol || !istr) {
    printf("Allocation failure for a\n");
    exit_status = -1;
    goto END;
  }

  /*  Concatenate sparse matrices B and C to form A */
  for (i = 0; i < bnnz; ++i) {
    a[i] = b[i];
    irow[i] = brow[i];
    icol[i] = bcol[i];
  }
  for (i = 0; i < cnnz; ++i) {
    a[bnnz + i] = c[i];
    irow[bnnz + i] = crow[i];
    icol[bnnz + i] = ccol[i];
  }

  /* Reorder along rows, sum duplicates and remove zeros */

  dup = Nag_SparseNsym_SumDups;
  zero = Nag_SparseNsym_RemoveZeros;

  /* nag_sparse_real_gen_sort (f11zac).
   * Sparse sort (nonsymmetric)
   */
  nag_sparse_real_gen_sort(n, &nnz, a, irow, icol, dup, zero, istr, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparse_real_gen_sort (f11zac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Output results */
  printf("Summed and reordered elements of A = B + C, along rows first\n");
  printf("nnz = %4" NAG_IFMT "\n", nnz);
  printf("%8s%16s%8s%8s\n", "", "a", "irow", "icol");

  for (i = 0; i < nnz; ++i)
    printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "%8" NAG_IFMT "\n", i, a[i],
           irow[i], icol[i]);

  /* Reorder down columns, fail on duplicates or zeros.
   * Creates CCS storage format as side-effect
   */
  dup = Nag_SparseNsym_FailDups;
  zero = Nag_SparseNsym_FailZeros;
  INIT_FAIL(fail);

  /* nag_sparse_real_gen_sort (f11zac).
   * Sparse sort (nonsymmetric)
   */
  nag_sparse_real_gen_sort(n, &nnz, a, icol, irow, dup, zero, istr, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparse_real_gen_sort (f11zac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Output results */
  printf("Reordered elements, along columns first\n");
  printf("nnz = %4" NAG_IFMT "\n", nnz);
  printf("%8s%16s%8s%8s\n", "", "a", "irow", "icol");

  for (i = 0; i < nnz; ++i)
    printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "%8" NAG_IFMT "\n", i, a[i],
           irow[i], icol[i]);

  /* Output results */
  printf("Same matrix in CCS format\n");
  printf("nnz = %4" NAG_IFMT "\n", nnz);
  printf("%8s%16s%8s%8s\n", "", "a", "irowix", "icolzp");

  for (i = 0; i < nnz; ++i) {
    if (i < n + 1) {
      printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "%8" NAG_IFMT "\n", i, a[i],
             irow[i], istr[i]);
    } else {
      printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "\n", i, a[i], irow[i]);
    }
  }

  /* Calculate 1-norm from Compressed Column Storage format */
  norm = Nag_RealOneNorm;
  INIT_FAIL(fail);

  /* nag_sparse_direct_real_gen_norm (f11mlc).
   * 1-norm, infinity-norm, largest absolute element, real
   * general matrix
   */
  nag_sparse_direct_real_gen_norm(norm, &anorm, n, istr, irow, a, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparse_direct_real_gen_norm (f11mlc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  /* Output norm */
  printf("%s %16.4e\n", "One-norm", anorm);

  /* Convert CCS format back to Coordinate Storage */
  icol2 = NAG_ALLOC(nnz, Integer);

  for (i = 0; i < n; ++i)
    for (j = istr[i] - 1; j < istr[i + 1] - 1; ++j)
      icol2[j] = i + 1;

  /* Check that result matches original Coordinate Storage Data */
  for (i = 0; i < nnz; ++i)
    if (icol[i] != icol2[i]) {
      printf("Converted CCS format back to CS format, result differs "
             "from original!\n");
      goto END;
    }
  printf("Converted CCS format back to CS format, result matches original\n");

  /* Swap the row and column vectors to obtain transpose of A */
  trrow = icol;
  trcol = irow;

  /* Output transposed results */
  printf("Transpose of summed and reordered elements, along rows first\n");
  printf("nnz = %4" NAG_IFMT "\n", nnz);
  printf("%8s%16s%8s%8s\n", "", "a", "irow", "icol");

  for (i = 0; i < nnz; ++i)
    printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "%8" NAG_IFMT "\n", i, a[i],
           trrow[i], trcol[i]);

END:
  NAG_FREE(istr);
  NAG_FREE(brow);
  NAG_FREE(c);
  NAG_FREE(crow);
  NAG_FREE(irow);
  NAG_FREE(bcol);
  NAG_FREE(ccol);
  NAG_FREE(icol);
  NAG_FREE(a);
  NAG_FREE(icol2);
  NAG_FREE(b);

  return exit_status;
}