subroutine matvec(job,a,ia,ja,b,x,y,n)
c****************************************************************************
c sparse matrix multiplication
c
c multiplies ==> ax = y
c   a = sparse matrix given by a,ia (and ja, if used)
c   x = vector given by x
c   y = vector returned by this routine with y = a*x
c   n = dimensions of a,x,y
c
c two possible storage schemes:
c
c scheme 1  -  matrix stored in a, ia, and ja:
c --------------------------------------------
c n   = number of variables/equations
c a   = contains the nonzero elements of a, stored by rows.
c       size = number of nonzero entries in a.
c ia  = pointers to delimit the rows in a.
c       size = n+1.
c       the index in ja and a of the first location following the last
c       element in the last row is stored in ia(n+1).
c ja  = column numbers corresponding to the elements of a.
c       size = size of a.
c
c this storage scheme is as follows:
c (1) ia(1),...,ia(n+1)
c         thus, ia(1),...,ia(n) will contain the positions in a in which the
c         first nonzero elements of rows 1,...,n are in a.
c (2) ia(n+1) contains the lengths of a and ja.
c         (points to first element past last in each)
c (3) ja(1),...,ja(ia(n+1)).
c         thus, ja(1),...,ja(ia(n+1)) will contain the column numbers
c         corresponding to the elements of a. (so size ja = size a).
c (4) a(1),...,a(ia(n+1)) are the nonzero elements of a.
c
c matrix-vector multiplication psuedocode interms of ia and ja:
c     *** do the multiplication ***
c     do 10 i = 1, n
c        y(i) = 0.0d0
c        do 10 j = ia(i), ia(i+1) - 1
c           y(i) = y(i) + a(j) * x(ja(j))
c10   continue
c
c scheme 2  -  matrix stored in a and iax:
c ---------------------------------------------
c n   = number of variables/equations
c a   = contains the nonzero elements of a, stored by rows.
c       size = number of nonzero entries in a.
c iax = (ia,ja), where ia and ja correspond to ia and ja in the sparse
c       matrix storage scheme used by ysmp.
c       ia = pointers to delimit the rows in a.
c          size = n+1.
c          the index in ja and a of the first location following the last
c          element in the last row is stored in ia(n+1).
c       ja = column numbers corresponding to the elements of a.
c          size = size of a.
c
c this storage scheme is as follows:
c (1) iax(1),...,iax(n+1) = ia(1),...,ia(n+1).
c         thus, ia(1),...,ia(n) will contain the positions in a in which the
c         first nonzero elements of rows 1,...,n are in a.
c (2) ia(n+1) contains the lengths of a and ja.
c         (points to first element past last in each)
c (3) iax(n+1+1),...,iax(n+1+iax(n+1)) = ja(1),...,ja(iax(n+1)).
c         thus, ja(1),...,ja(iax(n+1)) will contain the column numbers
c         corresponding to the elements of a. (so size ja = size a).
c (4) a(1),...,a(iax(n+1)) are the nonzero elements of a.
c
c matrix-vector multiplication psuedocode for iax, decoded into ia and ja:
c     *** decode iax into ia and ja ***
c     do 10 i = 1, n+1
c        ia(i) = iax(i)
c10   continue
c     length = ia(n+1) - 1
c     do 20 i = 1, length
c        ja(i) = iax(n+1+i)
c20   continue
c     *** do the multiplication ***
c     do 30 i = 1, n
c        y(i) = 0.0d0
c        do 30 j = ia(i), ia(i+1) - 1
c           y(i) = y(i) + a(j) * x(ja(j))
c30   continue
c
c matrix-vector multiplication psuedocode interms of iax directly:
c     *** alternate form ***
c     offset = n+1
c     do 10 i = 1, n
c        y(i) = 0.0d0
c        do 10 j = iax(i), iax(i+1) - 1
c           y(i) = y(i) + a(j) * x(iax(offset+j))
c10   continue
c
c author -->  michael jay holst
c date   -->  04 may 1989
c****************************************************************************
      complex a(*),b(*),x(*),y(*)
      integer ia(*),ja(*)
c
c     *** do the multiplication ***
      do 10 i = 1, n
         y(i) = (0.0,0.0)
         do 10 j = ia(i), ia(i+1) - 1
            y(i) = y(i) + a(j) * x(ja(j))
 10   continue
c
c     *** go home ***
      return
      end
c
      subroutine prtmat(a,ia,ja,n)
c****************************************************************************
      implicit  double precision(a-h,o-z)
      complex a(*)
      integer ia(*),ja(*)
c
c     *** do the printing ***
      do 10 i = 1, n
         do 10 j = ia(i), ia(i+1) - 1
            print*,' a(',i,',',ja(j),') = ',a(j)
 10   continue
c     *** go home ***
      return
      end