*=======================================================================
* program main
* test program for sparse operations
      program main
      implicit none
      integer m, n, nmbmax
      parameter (m=5, n=4, nmbmax=m*n)
      integer a_nmbelems, b_nmbelems, a_rowbase(0:m), b_rowbase(0:n),
     &  a_colind(0:nmbmax-1), b_colind(0:nmbmax-1)
      double precision a_elems(0:nmbmax-1), b_elems(0:nmbmax-1)
      double precision a_nrm1, a_nrmi, b_nrm1, b_nrmi
      double precision sp_nrm1, sp_nrmi
      external sp_setup, sp_transp, sp_nrm1, sp_nrmi
*     ------------------------------------------------------------------
*     ..define test problem
      call sp_setup (m, n, nmbmax, a_nmbelems, a_rowbase, a_colind,
     &  a_elems)
      write (*,*) 'a_nmbelems', a_nmbelems
      write (*,*) 'a_rowbase', a_rowbase
      write (*,*) 'a_colind', a_colind
      write (*,*) 'a_elems', a_elems
*     ..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)
*     ..form transpose
      b_nmbelems = a_nmbelems
      call sp_transp (m, n, a_nmbelems, a_rowbase, a_colind, a_elems,
     &  b_rowbase, b_colind, b_elems)
*     ..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)
*     ..write norms
      write (*,'(2x,a/4(4x,a,1pe12.2/))')
     &  'Matrix norms:',
     &  'a_nrm1 =', a_nrm1,
     &  'a_nrmi =', a_nrmi,
     &  'b_nrm1 =', b_nrm1,
     &  'b_nrmi =', b_nrmi
      stop
      end
*=======================================================================
* subroutine sp_setup
* define test matrix
      subroutine sp_setup (m, n, nmbmax, nmbelems, rowbase, colind,
     &  elems)
      implicit none
      integer m, n, nmbmax, nmbelems, rowbase(0:m), colind(0:nmbmax-1)
      double precision elems(0:nmbmax-1)
      integer ord, i, k
*     ------------------------------------------------------------------
*     The test matrix is an upper bi-diagonal matrix of order min(m,n).
*     For 0.le.i.lt.ord, A(i,i)=i+1;
*     For 0.le.i.lt.ord-1, A(i,i+1)=-1.
*     All other entries are zero.
      ord = min(m,n)
      if (nmbmax.lt.2*ord-1) then
       stop 'sp_setup -- declared dimension is insufficient'
      endif
      k = 0
      i = 0
      do while (i.lt.ord)
       rowbase(i) = k
       colind(k) = i
       elems(k) = i+1
       k = k+1
       if (i.lt.ord-1) then
        colind(k) = i+1
        elems(k) = -1
        k = k+1
       endif
       i = i+1
      enddo
      do while (i.le.m)
       rowbase(i) = k
       i = i+1
      enddo
      nmbelems = k
      return
      end
*=======================================================================
* double precision function sp_nrm1
* compute sparse matrix 1-norm
      double precision function sp_nrm1 (m, n, nmbelems, rowbase,
     &  colind, elems)
      implicit none
      integer m, n, nmbelems, rowbase(0:m), colind(0:nmbelems-1)
      double precision elems(0:nmbelems-1)
*     ------------------------------------------------------------------
      stop
      end
*=======================================================================
* double precision function sp_nrmi
* compute sparse matrix inf-norm
      double precision function sp_nrmi (m, n, nmbelems, rowbase,
     &  colind, elems)
      implicit none
      integer m, n, nmbelems, rowbase(0:m), colind(0:nmbelems-1)
      double precision elems(0:nmbelems-1)
*     ------------------------------------------------------------------
      stop
      end
*=======================================================================
* subroutine sp_trans
* form sparse matrix transpose
      subroutine sp_transp (m, n, nmbelems, a_rowbase, a_colind,
     &  a_elems, b_rowbase, b_colind, b_elems)
*     ..input arguments
      integer m, n, nmbelems, a_rowbase(0:m), a_colind(0:nmbelems-1)
      double precision a_elems(0:nmbelems-1)
*     ..output arguments
      integer b_rowbase(0:n), b_colind(0:nmbelems-1)
      double precision b_elems(0:nmbelems-1)
*     ------------------------------------------------------------------
      stop
      end