c fishpk38 from portlib                                  12/30/83
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     *                                                               *
c     *                        f i s h p a k                          *
c     *                                                               *
c     *                                                               *
c     *     a package of fortran subprograms for the solution of      *
c     *                                                               *
c     *      separable elliptic partial differential equations        *
c     *                                                               *
c     *                  (version 3.1 , october 1980)                  *
c     *                                                               *
c     *                             by                                *
c     *                                                               *
c     *        john adams, paul swarztrauber and roland sweet         *
c     *                                                               *
c     *                             of                                *
c     *                                                               *
c     *         the national center for atmospheric research          *
c     *                                                               *
c     *                boulder, colorado  (80307)  u.s.a.             *
c     *                                                               *
c     *                   which is sponsored by                       *
c     *                                                               *
c     *              the national science foundation                  *
c     *                                                               *
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     program to illustrate the use of subroutine genbun to
c     solve the equation
c
c     (1+x)**2*(d/dx)(du/dx) - 2(1+x)(du/dx) + (d/dy)(du/dy)
c
c                  = 3(1+x)**4*sin(y)                      (1)
c
c     on the rectangle 0 .lt. x .lt. 1 and -pi .lt. y .lt. pi
c     with the boundary conditions
c
c     (du/dx)(0,y) = 4sin(y)                               (2)
c                                -pi .le. y .le. pi
c     u(1,y) = 16sin(y)                                    (3)
c
c     and with u periodic in y using finite differences on a
c     grid with deltax (= dx) = 1/20 and deltay (= dy) = pi/20.
c        to set up the finite difference equations we define
c     the grid points
c
c     x(i) = (i-1)dx            i=1,2,...,21
c
c     y(j) = -pi + (j-1)dy      j=1,2,...,41
c
c     and let v(i,j) be an approximation to u(x(i),y(j)).
c     numbering the grid points in this fashion gives the set
c     of unknowns as v(i,j) for i=1,2,...,20 and j=1,2,...,40.
c     hence, in the program m = 20 and n = 40.  at the interior
c     grid point (x(i),y(j)), we replace all derivatives in
c     equation (1) by second order central finite differences,
c     multiply by dy**2, and collect coefficients of v(i,j) to
c     get the finite difference equation
c
c     a(i)v(i-1,j) + b(i)v(i,j) + c(i)v(i+1,j)
c
c     + v(i,j-1) - 2v(i,j) + v(i,j+1) = f(i,j)            (4)
c
c     where s = (dy/dx)**2, and for i=2,3,...,19
c
c     a(i) = (1+x(i))**2*s + (1+x(i))*s*dx
c
c     b(i) = -2(1+x(i))**2*s
c
c     c(i) = (1+x(i))**2*s - (1+x(i))*s*dx
c
c     f(i,j) = 3(1+x(i))**4*dy**2*sin(y(j))  for j=1,2,...,40.
c
c        to obtain equations for i = 1, we replace the
c     derivative in equation (2) by a second order central
c     finite difference approximation, use this equation to
c     eliminate the virtual unknown v(0,j) in equation (4)
c     and arrive at the equation
c
c     b(1)v(1,j) + c(1)v(2,j) + v(1,j-1) - 2v(1,j) + v(1,j+1)
c
c                       = f(1,j)
c
c     where
c
c     b(1) = -2s , c(1) = 2s
c
c     f(1,j) = (11+8/dx)*dy**2*sin(y(j)),  j=1,2,...,40.
c
c     for completeness, we set a(1) = 0.
c        to obtain equations for i = 20, we incorporate
c     equation (3) into equation (4) by setting
c
c     v(21,j) = 16sin(y(j))
c
c     and arrive at the equation
c
c     a(20)v(19,j) + b(20)v(20,j)
c
c     + v(20,j-1) - 2v(20,j) + v(20,j+1) = f(20,j)
c
c     where
c
c     a(20) = (1+x(20))**2*s + (1+x(20))*s*dx
c
c     b(20) = -2*(1+x(20))**2*s
c
c     f(20,j) = (3(1+x(20))**4*dy**2 - 16(1+x(20))**2*s
c                + 16(1+x(20))*s*dx)*sin(y(j))
c
c                    for j=1,2,...,40.
c
c     for completeness, we set c(20) = 0.  hence, in the
c     program mperod = 1.
c        the periodicity condition on u gives the conditions
c
c     v(i,0) = v(i,40) and v(i,41) = v(i,1) for i=1,2,...,20.
c
c     hence, in the program nperod = 0.
c
c          the exact solution to this problem is
c
c                  u(x,y) = (1+x)**4*sin(y) .
c
      dimension       f(22,40)   ,a(20)      ,b(20)      ,c(20)      ,
     1                x(21)      ,y(41)      ,w(380)
c
c     from the dimension statement we get that idimf = 22 and that w
c     has been dimensioned
c
c     4n + (10+int(log2(n)))m = 4*20 + (10+5)*20 = 380 .
c
      idimf = 22
      m = 20
      mp1 = m+1
      mperod = 1
      dx = 0.05
      n = 40
      nperod = 0
      pi = pimach(dum)
      dy = pi/20.
c
c     generate grid points for later use.
c
      do 101 i=1,mp1
         x(i) = float(i-1)*dx
  101 continue
      do 102 j=1,n
         y(j) = -pi+float(j-1)*dy
  102 continue
c
c     generate coefficients.
c
      s = (dy/dx)**2
      do 103 i=2,19
         t = 1.+x(i)
         tsq = t**2
         a(i) = (tsq+t*dx)*s
         b(i) = -2.*tsq*s
         c(i) = (tsq-t*dx)*s
  103 continue
      a(1) = 0.
      b(1) = -2.*s
      c(1) = -b(1)
      b(20) = -2.*s*(1.+x(20))**2
      a(20) = -b(20)/2.+(1.+x(20))*dx*s
      c(20) = 0.
c
c     generate right side.
c
      do 105 i=2,19
         do 104 j=1,n
            f(i,j) = 3.*(1.+x(i))**4*dy**2*sin(y(j))
  104    continue
  105 continue
      t = 1.+x(20)
      tsq = t**2
      t4 = tsq**2
      do 106 j=1,n
         f(1,j) = (11.+8./dx)*dy**2*sin(y(j))
         f(20,j) = (3.*t4*dy**2-16.*tsq*s+16.*t*s*dx)*sin(y(j))
  106 continue
      call genbun (nperod,n,mperod,m,a,b,c,idimf,f,ierror,w)
c
c     compute discretiazation error.  the exact solution is
c
c            u(x,y) = (1+x)**4*sin(y) .
c
      err = 0.
      do 108 i=1,m
         do 107 j=1,n
            t = abs(f(i,j)-(1.+x(i))**4*sin(y(j)))
            if (t .gt. err) err = t
  107    continue
  108 continue
      print 1001 , ierror,err,w(1)
      stop
c
 1001 format (1h1,20x,25hsubroutine genbun example///
     1        10x,46hthe output from the ncar control data 7600 was//
     2        32x,10hierror = 0/
     3        18x,34hdiscretization error = 9.64063e-03/
     4        12x,32hrequired length of w array = 380//
     5        10x,32hthe output from your computer is//
     6        32x,8hierror =,i2/18x,22hdiscretization error =,e12.5/
     7        12x,28hrequired length of w array =,f4.0)
c
      end
c     compute discretiazation error.  the exact solution is