#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

struct spmat {
  int m, n, nmbelems, *rowbase, *colind;
  double *elems;
};

/* sp_nrmi: return sparse matrix inf-norm.
   sp_nrmi(aa) = (Max i: 0<=i<m: (Sum j: 0<=j<n: fabs(A(i,j)))). */
double sp_nrmi(struct spmat aa)
{
  double t0, t1; int i, k;
  t0 = 0; i = 0;
  /* loop invariant:
       t0 = (Max ii: 0<=ii<i: (Sum j: 0<=j<n: fabs(A(i,j)))).
     Because we are dealing with non-negative values only we take
     the max of the empty set as 0. */
  while (i!=aa.m) {
    /* compute row-sum for row i */
    t1 = 0; k = aa.rowbase[i];
    while (k!=aa.rowbase[i+1]) {
      t1 = t1+fabs(aa.elems[k]);
      ++k;
    }
    if (t0<t1) {
      t0 = t1;
    }
    ++i;
  }
  /* (loop invariant) and (i==aa.m) hence
     t0 = (Max ii: 0<=ii<m: (Sum j: 0<=j<n: fabs(A(i,j)))) */
  return t0;
}

/* sp_nrm1: return sparse matrix 1-norm.
   sp_nrm1(aa) = (Max j: 0<=j<n: (Sum i: 0<=i<m: fabs(A(i,j)))). */
double sp_nrm1(struct spmat aa)
{
  double *colnrm; double t0; int j, k;
  /* colnrm is intended as [0:n-1] double array; it will accumulate the
     column-sums */
  colnrm = (double *) malloc(aa.n*sizeof(double));
  for (j=0; j<aa.n; j++) {
    colnrm[j] = 0;
  }
  for (k=0; k<aa.nmbelems; k++) {
    colnrm[aa.colind[k]] = colnrm[aa.colind[k]]+fabs(aa.elems[k]);
  }
  /* compute max elem in colnrm */
  t0 = 0; j = 0;
  for (j=0; j<aa.n; j++) {
    if (t0<colnrm[j]) {
      t0 = colnrm[j];
    }
  }
  return t0;
}

/* sp_transp: form sparse matrix transpose. */
struct spmat sp_transp(struct spmat aa)
{
  struct spmat bb;
  int *colcount; int i, j, k, k1;
  bb.m = aa.n; bb.n = aa.m; bb.nmbelems = aa.nmbelems;
  /* storage allocation for colcount: [0:n-1] int workspace. */
  colcount = (int *) malloc(aa.n*sizeof(int));
  /* storage allocation for arrays in bb */
  bb.rowbase = (int *) malloc((aa.n+1)*sizeof(int));
  bb.colind = (int *) malloc(aa.nmbelems*sizeof(int));
  bb.elems = (double *) malloc(aa.nmbelems*sizeof(double));
  /* count number of elements in each column of aa */
  for (j=0; j<aa.n; j++) {
    colcount[j] = 0;
  }
  for (k=0; k<aa.nmbelems; k++) {
    colcount[aa.colind[k]] = colcount[aa.colind[k]]+1;
  }
  /* compute bb.rowbase */
  bb.rowbase[0] = 0;
  for (j=0; j<aa.n; j++) {
    bb.rowbase[j+1] = bb.rowbase[j]+colcount[j];
  }
  /* we next assign bb.colind and bb.elems, running through the
     elements in the natural order for aa (which is a scattered order
     for bb).  For 0<=j<n, colcount[j] will accumulate the number of
     elements visited in column j of aa. */
  for (j=0; j<aa.n; j++) {
    colcount[j] = 0;
  }
  for (i=0; i<aa.m; i++) {
    for (k=aa.rowbase[i]; k!=aa.rowbase[i+1]; k++) {
      j = aa.colind[k];
      /* compute the index k1 for aa's k-th element in the array bb */
      k1 = bb.rowbase[j]+colcount[j];
      /* assign bb.colind and bb.elems for this element */
      bb.colind[k1] = i;
      bb.elems[k1] = aa.elems[k];
      /* update colcount */
      colcount[j] = colcount[j]+1;
    }
  }
  return bb;
}

/* sp_setup: define sparse test matrix. */
struct spmat sp_setup(int m, int n)
{
  struct spmat aa;
  /* The test matrix is an upper bi-diagonal matrix of order min(m,n).
     For 0<=i<ord, A(i,i)=i+1;
     For 0<=i<ord-1, A(i,i+1)=-1. */
  int ord, k, i;
  if (m<n) {
    ord = m;
  }
  else {
    ord = n;
  }
  aa.m = m; aa.n = n;
  aa.nmbelems = 2*ord-1;
  aa.rowbase = (int *) malloc((aa.m+1)*sizeof(int));
  aa.colind = (int *) malloc(aa.nmbelems*sizeof(int));
  aa.elems = (double *) malloc(aa.nmbelems*sizeof(double));
  k = 0; i = 0;
  while (i<ord) {
    aa.rowbase[i] = k;
    aa.colind[k] = i;
    aa.elems[k] = i+1;
    ++k;
    if (i<ord-1) {
      aa.colind[k] = i+1;
      aa.elems[k] = -1;
      ++k;
    }
    ++i;
  }
  while (i<=m) {
    aa.rowbase[i] = k;
    ++i;
  }
  assert(k==aa.nmbelems);
  return aa;
}

/* sp_setup1: define sparse test matrix. */
struct spmat sp_setup1(int m, int n)
{
  struct spmat aa;
  /* The test matrix is an upper bi-diagonal matrix of order min(m,n)
     with one additional element in the lower left corner.
     For 0<=i<ord, A(i,i)=i+1;
     For 0<=i<ord-1, A(i,i+1)=-1;
     A(ord-1,0)=-0.5 */
  int ord, k, i;
  if (m<n) {
    ord = m;
  }
  else {
    ord = n;
  }
  aa.m = m; aa.n = n;
  aa.nmbelems = 2*ord;
  aa.rowbase = (int *) malloc((aa.m+1)*sizeof(int));
  aa.colind = (int *) malloc(aa.nmbelems*sizeof(int));
  aa.elems = (double *) malloc(aa.nmbelems*sizeof(double));
  k = 0; i = 0;
  while (i<ord) {
    aa.rowbase[i] = k;
    if (i==ord-1) {
      aa.colind[k] = 0;
      aa.elems[k] = -0.5;
      ++k;
    }
    aa.colind[k] = i;
    aa.elems[k] = i+1;
    ++k;
    if (i<ord-1) {
      aa.colind[k] = i+1;
      aa.elems[k] = -1;
      ++k;
    }
    ++i;
  }
  while (i<=m) {
    aa.rowbase[i] = k;
    ++i;
  }
  assert(k==aa.nmbelems);
  return aa;
}
 
main()
{
  int m, n, i, k;
  struct spmat aa, bb;
  /* form test matrix aa */
  m = 5; n = 4;
  aa = sp_setup1 (m, n);
  /* print aa */
  printf("  %s %4d  %s %4d  %s %4d\n",
	 "aa.m", aa.m, "aa.n", aa.n, "aa.nmbelems", aa.nmbelems);
  printf("  aa.rowbase\n");
  for (i=0; i<=aa.m; i++) {
    printf(" %4d", aa.rowbase[i]);
  }
  printf("\n  aa.colind\n");
  for (k=0; k<aa.nmbelems; k++) {
    printf(" %4d", aa.colind[k]);
  }
  printf("\n  aa.elems\n");
  for (k=0; k<aa.nmbelems; k++) {
    printf(" %12.2e", aa.elems[k]);
  }
  printf("\n");
  /* print norms of aa */
  printf("  %s  %12.2e    %s  %12.2e\n",
    "aa.nrm1", sp_nrm1(aa), "aa.nrmi", sp_nrmi(aa));
  /* form transpose matrix */
  bb = sp_transp(aa);
  /* print bb */
  printf("  %s %4d  %s %4d  %s %4d\n",
	 "bb.m", bb.m, "bb.n", bb.n, "bb.nmbelems", bb.nmbelems);
  printf("  bb.rowbase\n");
  for (i=0; i<=bb.m; i++) {
    printf(" %4d", bb.rowbase[i]);
  }
  printf("\n  bb.colind\n");
  for (k=0; k<bb.nmbelems; k++) {
    printf(" %4d", bb.colind[k]);
  }
  printf("\n  bb.elems\n");
  for (k=0; k<bb.nmbelems; k++) {
    printf(" %12.2e", bb.elems[k]);
  }
  printf("\n");
  /* print norms of bb */
  printf("  %s  %12.2e    %s  %12.2e\n",
    "bb.nrm1", sp_nrm1(bb), "bb.nrmi", sp_nrmi(bb));
}