program main
c
c ... array declarations.
c
      real    rhs(100), u(100), wksp(600), ubar(1), rparm(30),
     a        coef(1), wfac(1)
      integer iparm(30), jcoef(1), jwfac(1)
      external mult, srpass
c
      nw = 600
c
c ... generate rhs. 
c
      nx = 10
      ny = 10
      n = nx*ny
      h = 1.0/float(nx + 1)
      do 10 i = 1,n 
         rhs(i) = 0.0
 10   continue
      k = 0
      do 30 j = 1,ny
         y = float(j)*h
         do 25 i = 1,nx
            x = float(i)*h
            k = k + 1
            if (j .eq. 1) rhs(k) = rhs(k) + 2.0
            if (j .eq. ny) rhs(k) = rhs(k) + 2.0*(1.0 + x)
            if (i .eq. 1) rhs(k) = rhs(k) + 1.0
            if (i .eq. nx) rhs(k) = rhs(k) + 1.0 + y
 25      continue
 30   continue
      call dfault (iparm,rparm)
c
c ... now, reset some default values.
c
      iparm(3) = 3
      rparm(1) = 1.0e-8
c
c ... generate an initial guess for u and call nspcg.
c
      call vfill (n,u,0.0)
c
      call sorw (mult,srpass,coef,jcoef,wfac,jwfac,n,
     a            u,ubar,rhs,wksp,nw,iparm,rparm,ier)
      stop
      end 
      subroutine mult (coef,jcoef,wfac,jwfac,n,x,y)
      real x(n), y(n)
      nx = 10
c
c ... compute product as if first superdiagonal and first
c     subdiagonal were full.
c
      y(1) = 6.0*x(1) - x(2) - 2.0*x(nx+1)
      do 10 i = 2,nx
         y(i) = 6.0*x(i) - x(i+1) - x(i-1) - 2.0*x(i+nx)
 10   continue
      do 15 i = nx+1,n-nx
         y(i) = 6.0*x(i) - x(i+1) - x(i-1) - 2.0*(x(i+nx) + x(i-nx))
 15   continue
      do 20 i = n-nx+1,n-1
         y(i) = 6.0*x(i) - x(i+1) - x(i-1) - 2.0*x(i-nx)
 20   continue
      y(n) = 6.0*x(n) - x(n-1) - 2.0*x(n-nx)
c
c ... make corrections to y vector for zeros in first superdiagonal
c     and first subdiagonal.
c
      do 25 i = nx,n-nx,nx
         y(i) = y(i) + x(i+1) 
 25   continue
      do 30 i = nx+1,n-nx+1,nx
         y(i) = y(i) + x(i-1) 
 30   continue
      return
      end 
      subroutine srpass (coef,jcoef,wfac,jwfac,n,u,rhs,unew)
c
c ... srpass does an sor iteration.
c
c        unew = inv((1/w)*d + l)*(((1/w)-1)*d*un + rhs - u*un)
c
c ... parameters -- 
c
c        n      order of system
c        u      current solution vector 
c        rhs    right hand side
c        unew   updated solution vector 
c
c ... specifications for parameters
c
      dimension coef(1), wfac(1), u(1), unew(1), rhs(1)
      integer   jcoef(1), jwfac(1)
      logical omgadp
      common / itcom5 / omgadp, omega, alphab, betab, fff, specr
c
c ... temp = ((1/w)-1)*d*un + rhs - u*un
c     unew is used for temp.
c
      nx = 10
      con = 6.0*(1.0/omega - 1.0)
      do 10 i = 1,n 
         unew(i) = con*u(i) + rhs(i)
 10   continue
      do 15 i = 1,n-1
         unew(i) = unew(i) + u(i+1)
 15   continue
      do 20 i = 1,n-nx
         unew(i) = unew(i) + 2.0*u(i+nx)
 20   continue
      do 25 i = nx,n-nx,nx
         unew(i) = unew(i) - u(i+1)
 25   continue
c
c ... unew = inv((1/w)*d + l)*temp
c
      con = omega/6.0
      do 40 j = 1,nx
         ibgn = (j-1)*nx + 1
         iend = j*nx
         unew(ibgn) = con*unew(ibgn)
         do 30 i = ibgn+1,iend
            unew(i) = con*(unew(i) + unew(i-1))
 30      continue
         if (j .eq. nx) go to 40
cdir$ ivdep
         do 35 i = ibgn,iend
            unew(i+nx) = unew(i+nx) + 2.0*unew(i) 
 35      continue
 40   continue
      return
      end