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

NAG CL Interface Introduction
Example description
/* nag_sparse_real_rect_sort (f11zcc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.3, 2022.
 */
/* Pre-processor includes */
#include <nag.h>
#include <stdio.h>

int main(void) {
  /* Double scalar and array declarations */
  double tmpa;
  double *a = 0, *b = 0, *c = 0;

  /* Integer scalar and array declarations. */
  Integer bnnz, cnnz, i, m, n, nnz, tmprow;
  Integer exit_status = 0;
  Integer *bcol = 0, *brow = 0, *ccol = 0, *crow = 0;
  Integer *icol = 0, *irow = 0, *istc = 0;

  /* NAG structures and types */
  Nag_SparseNsym_Dups dup;
  Nag_SparseNsym_Store store;
  Nag_SparseNsym_Zeros zer;
  NagError fail;

  /* Initialize the error structure. */
  INIT_FAIL(fail);

  printf("nag_sparse_real_rect_sort (f11zcc) Example Program Results\n\n");

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

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

  /* Allocate memory for B */
  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 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 entries for C */
  scanf("%" NAG_IFMT "%*[^\n] ", &cnnz);

  /* Allocate memory for C */
  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 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 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("\n");

  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]);
  }
  printf("\n");

  /* Set initial size for A, and allocate memory */
  nnz = bnnz + cnnz;
  a = NAG_ALLOC(nnz, double);
  irow = NAG_ALLOC(nnz, Integer);
  icol = NAG_ALLOC(nnz, Integer);
  istc = NAG_ALLOC(n + 1, Integer);

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

  /* Concatentate 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];
  }

  /* A is in Coordinate Storage (CS) format */
  store = Nag_SparseNsym_StoreCS;

  /* Sum duplicates and remove zeros */
  dup = Nag_SparseNsym_SumDups;
  zer = Nag_SparseNsym_RemoveZeros;

  /* Sort into column major order */
  nag_sparse_real_rect_sort(m, n, &nnz, a, irow, icol, istc, store, dup, zer,
                            &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_spare_real_rect_sort (f11zcc)\n");
    printf("%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Output results */
  printf("Summed and reordered elements of A = B + C, ");
  printf("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]);
  }
  printf("\n");

  /* Output results in Compressed Column Storage (CCS) format */
  printf("Same matrix in CCS format\n");
  printf("NNZ = %4" NAG_IFMT "\n", nnz);
  printf("%8s%16s%8s%8s\n", "", "A", "IROW", "ISTC");
  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], istc[i]);
    } else {
      printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "\n", i, a[i], irow[i]);
    }
  }
  printf("\n");

  /* Reorder some rows within columns */

  /* 1,1 <--> 3,1 at i=0,1 */
  tmpa = a[0];
  tmprow = irow[0];
  a[0] = a[1];
  irow[0] = irow[1];
  a[1] = tmpa;
  irow[1] = tmprow;

  /* 1,3 <--> 3,3 at i=4,6 */
  tmpa = a[4];
  tmprow = irow[4];
  a[4] = a[6];
  irow[4] = irow[6];
  a[6] = tmpa;
  irow[6] = tmprow;

  /* Output (in CCS format) the matrix being converted */
  printf("Same matrix with some rows put out of order\n");
  printf("NNZ = %4" NAG_IFMT "\n", nnz);
  printf("%8s%16s%8s%8s\n", "", "A", "IROW", "ISTC");
  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], istc[i]);
    } else {
      printf("%8" NAG_IFMT "%16.4e%8" NAG_IFMT "\n", i, a[i], irow[i]);
    }
  }
  printf("\n");

  /* Set up for conversion back to CS */
  for (i = 0; i < nnz; i++) {
    icol[i] = 0;
  }
  store = Nag_SparseNsym_StoreCCS;
  dup = Nag_SparseNsym_FailDups;
  zer = Nag_SparseNsym_FailZeros;
  INIT_FAIL(fail);

  nag_sparse_real_rect_sort(m, n, &nnz, a, irow, icol, istc, store, dup, zer,
                            &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_spare_real_rect_sort (f11zcc)\n");
    printf("%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Output in CS to compare to original */
  printf("Rows reordered and matrix converted from CCS to CS format\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]);
  }
  printf("\n");

  /* Use CS format to pass transpose of A */
  for (i = 0; i <= n; i++) {
    istc[i] = 0;
  }
  store = Nag_SparseNsym_StoreCS;
  INIT_FAIL(fail);

  /* Pass icol and irow in opposite order */
  nag_sparse_real_rect_sort(n, m, &nnz, a, icol, irow, istc, store, dup, zer,
                            &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_spare_real_rect_sort (f11zcc)\n");
    printf("%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Output in result in CS format */
  printf("Transpose of this matrix passed to get row-major order\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]);
  }
  printf("\n");

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

  return (exit_status);
}