#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)
{
  abort;
}

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