function spmat
% test program for sparse matrix manipulation.
  m = 5; n = 4;
  % define test matrix
  [a_nmbelems, a_rowbase, a_colind, a_elems] = sp_setup1 (m, n);
  % print a_rowbase and a_colind
  a_rowbase
  a_colind
  % compute norms
  a_nrm1 = sp_nrm1 (m, n, a_nmbelems, a_rowbase, a_colind, a_elems);
  a_nrmi = sp_nrmi (m, n, a_nmbelems, a_rowbase, a_colind, a_elems);
  % print norms of a
  fprintf('  %s  %12.2e    %s  %12.2e\n', 'nrm1', a_nrm1, 'nrmi', a_nrmi);
  % form transpose
  b_nmbelems = a_nmbelems;
  [b_rowbase, b_colind, b_elems] = sp_transp (m, n, a_nmbelems,...
    a_rowbase, a_colind, a_elems);
  % print b_rowbase and b_colind
  b_rowbase
  b_colind
  % compute norms of transpose
  b_nrm1 = sp_nrm1 (n, m, b_nmbelems, b_rowbase, b_colind, b_elems);
  b_nrmi = sp_nrmi (n, m, b_nmbelems, b_rowbase, b_colind, b_elems);
  % print norms of b
  fprintf('  %s  %12.2e    %s  %12.2e\n', 'nrm1', b_nrm1, 'nrmi', b_nrmi);

  function nrmi = sp_nrmi (m, n, nmbelems, rowbase, colind, elems)
% return sparse matrix inf-norm.
% nrmi = (Max i: 0<=i<m: (Sum j: 0<=j<n: abs(A(i,j)))).
  t0 = 0;
  % loop invariant:
  %   t0 = (Max ii: 1<=ii<i: (Sum j: 1<=j<=n: abs(A(i,j)))).
  % Because we are dealing with non-negative values only we take
  % the max of the empty set as 0.
  for i = 1:m
    % compute row-sum for row i.
    t1 = 0;
    for k = rowbase(i):rowbase(i+1)-1
      t1 = t1+abs(elems(k));
    end
    t0 = max(t0,t1);
  end
  nrmi = t0;

  function nrm1 = sp_nrm1 (m, n, nmbelems, rowbase, colind, elems)
% return sparse matrix 1-norm.
% nrm1 = (Max j: 0<=j<n: (Sum i: 0<=i<m: abs(A(i,j)))).
  colnrm = zeros(1,n);
  % colnrm will accumulate the column sums
  for k = 1:nmbelems
   colnrm(colind(k)) = colnrm(colind(k))+abs(elems(k));
  end
  nrm1 = max(colnrm);

  function [b_rowbase, b_colind, b_elems] = sp_transp (m, n, nmbelems,...
    rowbase, colind, elems)
% form transpose of sparse matrix
  error;

  function [nmbelems, rowbase, colind, elems] = sp_setup (m, n)
% set up sparse test matrix
% The test matrix is an upper bi-diagonal matrix of order min(m,n).
% For 1<=i<=ord, A(i,i)=i;
% For 1<=i<=ord-1, A(i,i+1)=-1.
  ord = min(m,n);
  nmbelems = 2*ord-1;
  rowbase = zeros(1,m+1);
  colind = zeros(1,nmbelems);
  elems = zeros(1,nmbelems);
  k = 1;
  for i = 1:ord
   rowbase(i) = k;
   colind(k) = i;
   elems(k) = i;
   k = k+1;
   if (i<ord)
    colind(k) = i+1;
    elems(k) = -1;
    k = k+1;
   end
  end
  for i = ord+1:m+1
   rowbase(i) = k;
  end
  if (k~=nmbelems+1)
    error ('wrong number of elements')
  end

  function [nmbelems, rowbase, colind, elems] = sp_setup1 (m, n)
% set up sparse test matrix
% The test matrix is an upper bi-diagonal matrix of order min(m,n)
% with one additional element in the lower left corner.
% For 1<=i<=ord, A(i,i)=i;
% For 1<=i<=ord-1, A(i,i+1)=-1;
% A(ord-1,0)=-0.5.
  ord = min(m,n);
  nmbelems = 2*ord;
  rowbase = zeros(1,m+1);
  colind = zeros(1,nmbelems);
  elems = zeros(1,nmbelems);
  k = 1;
  for i = 1:ord
   rowbase(i) = k;
   if (i==ord)
     colind(k) = 1;
     elems(k) = -0.5;
     k = k+1;
   end
   colind(k) = i;
   elems(k) = i;
   k = k+1;
   if (i<ord)
    colind(k) = i+1;
    elems(k) = -1;
    k = k+1;
   end
  end
  for i = ord+1:m+1
   rowbase(i) = k;
  end
  if (k~=nmbelems+1)
    error ('wrong number of elements')
  end