program main
c
c ... array declarations.
c
      real coef(120,4), rhs(100), u(100), wksp(600), ubar(1),
     a     rparm(30)
      integer jcoef(5), iwksp(300), iparm(30), p(1), ip(1)
      external cg, mic2
c
      ndim = 120
      mdim = 4
      n = 100
      maxnz = 3
      nw = 600
      inw = 300
c
      h = 1.0/11.0
c
c ... generate laplace*s equation.
c
      call matgen (n,ndim,maxnz,jcoef,coef,rhs,h,0)
      call dfault (iparm,rparm)
c
c ... now, reset some default values.
c
      iparm(2) = 50 
      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 nspcg (mic2,cg,ndim,mdim,n,maxnz,coef,jcoef,p,ip,
     a            u,ubar,rhs,wksp,iwksp,nw,inw,iparm,rparm,ier)
      stop
      end 
      subroutine matgen (n,ndim,maxnz,jcoef,coef,rhs,hh,isym)
c
c     matgen generates the itpack matrices jcoef, coef, and rhs
c     for a 5-point central difference approximation of the 
c     self-adjoint pde
c
c          (a*u )     +  (c*u )   + f*u  =  g
c              x x           y y
c
c     where the coefficients a, c, f, and g are functions of (x,y)
c     (supplied by the user) and the domain is the unit square
c     (0,1) x (0,1).  dirichlet boundary conditions are imposed
c     upon the boundary in the form of the function ub(x,y) also
c     supplied by the user.
c
c     parameters -- 
c
c        n       number of linear equations (output)
c        ndim    row dimension of jcoef and coef in the calling routine
c                 (input)
c        maxnz   maximum number of nonzeros per row (output)
c        jcoef   array of column indices (output) 
c        coef    array of equation coefficients (output)
c        rhs     vector of right-hand-side values (output)
c        h       mesh spacing (= hh) (input)
c
c     specifications for parameters
c
c     implicit half precision (a-h,o-z) 
c     implicit double precision (a-h,o-z)
      integer   jcoef(5)
      dimension coef(ndim,5), rhs(1)
      logical symm
c
c     statement functions
c
      a(x,y) = 1.0e0
      c(x,y) = 2.0e0
      f(x,y) = 0.0e0
      g(x,y) = 0.0e0
      ub(x,y) = 1.0e0 + x*y
c
      h = hh
      zero = 0.0e0
      one = 1.0e0
      two = 2.0e0
      tenth = 0.1e0 
      m = int ( one/h + tenth)
      mm1 = m - 1
      n = mm1*mm1
      halfh = h/two 
      symm   = isym .eq. 0
      maxnz = 5
      if (symm) maxnz = 3
c
c     loop on equations, assuming a natural ordering from left to right
c     and from down to up.
c
      neq = 0
      do 50 j = 1,mm1
         yy = float (j)*h
         do 40 i = 1,mm1
            xx = float (i)*h
            neq = neq + 1
            ae = a( xx+halfh, yy)
            aw = a( xx-halfh, yy)
            cn = c( xx, yy+halfh)
            cs = c( xx, yy-halfh)
            fp = f( xx, yy)
            gp = g( xx, yy)
            cc = ae + cn + aw + cs - h*h*fp
c
c ... center point
c
            coef(neq,1) = cc
            rhs(neq) = -h*h*gp
c
c ... east point
c
            if (i .eq. mm1) go to 5
            coef(neq,2) = -ae 
            go to 10
 5          coef(neq,2) = zero
            rhs(neq) = rhs(neq) + ae*ub( one, yy) 
c
c ... north point
c
 10         if (j .eq. mm1) go to 15
            coef(neq,3) = -cn 
            go to 20
 15         coef(neq,3) = zero
            rhs(neq) = rhs(neq) + cn*ub( xx, one) 
c
c ... west point
c
 20         if (i .eq. 1) go to 25
            if (.not. symm) coef(neq,4) = -aw
            go to 30
 25         if (.not. symm) coef(neq,4) = zero
            rhs(neq) = rhs(neq) + aw*ub(zero,yy)
c
c ... south point
c
 30         if (j .eq. 1) go to 35
            if (.not. symm) coef(neq,5) = -cs
            go to 40
 35         if (.not. symm) coef(neq,5) = zero
            rhs(neq) = rhs(neq) + cs*ub(xx,zero)
c
 40      continue
 50   continue
c
c ... data structure 2.
c
      jcoef(1) = 0
      jcoef(2) = 1
      jcoef(3) = mm1
      return
      end