!  fishprb2.f  23 February 1998
!
      program fishprb
!
!***********************************************************************
!
!! FISHPRB calls the FISHPACK tests.
!
      call test16

      stop
      end
      subroutine test16
!
!***********************************************************************
!
!! TEST16 uses BLKTRI to solve a simple PDE.
!
      integer m
      integer n

      parameter ( m = 11 )
      parameter ( n = 11 )
!
      integer idimy
      parameter ( idimy = m )
!
      real am(m)
      real an(n)
      real bm(m)
      real bn(n)
      real cm(m)
      real cn(n)
      real err
      real errmax
      integer i
      integer ierror
      integer iflg
      integer imax
      integer ipass
      integer j
      integer jmax
      integer mp
      integer np
      real w(1952)
      real x(m)
      real y(n)
      real yy(idimy,n)
!
      np = 1
      mp = 1

      delta = 0.1
!
!  Generate and store the grid points.
!
      do i = 1, m
        x(i) = real ( i-1 ) / real ( m-1 )
      end do

      do j = 1, n
        y(j) = real ( j-1) / real ( n-1 )
      end do
!
!  Compute the coefficients AM, BM, CM corresponding to the X direction.
!
      do i = 2, m-1
        am(i) = 1.0
        bm(i) = -2.0
        cm(i) = 1.0
      end do
!
!  Compute the coefficients AN, BN, CN corresponding to the Y direction.
!
      do j = 2, n-1
        an(j) = 1.0
        bn(j) = -2.0
        cn(j) = 1.0
      end do
!
!  Set the right hand side of the equations in the interior.
!
      do j = 2, n-1
        do i = 2, m-1
          yy(i,j) = 4.0 * delta**2
        end do
      end do
!
!  Set the boundary condition coefficients.
!
      am(1) = 0.0
      bm(1) = 1.0
      cm(1) = 0.0

      an(1) = 0.0
      bn(1) = 1.0
      cn(1) = 0.0

      am(m) = 0.0
      bm(m) = 1.0
      cm(m) = 0.0

      an(n) = 0.0
      bn(n) = 1.0
      cn(n) = 0.0
!
!  Assign the right hand side of boundary conditions.
!
      j = 1
      do i = 1, m
        yy(i,j) = x(i)**2 + y(j)**2 
      end do

      i = 1
      do j = 1, n
        yy(i,j) = x(i)**2 + y(j)**2
      end do

      i = m
      do j = 1, n
        yy(i,j) = x(i)**2 + y(j)**2
      end do

      j = n
      do i = 1, m
        yy(i,n) = x(i)**2 + y(j)**2
      end do

      yy(1,1) = 2.0 * yy(1,1)
      yy(m,1) = 2.0 * yy(m,1)
      yy(1,n) = 2.0 * yy(1,n)
      yy(m,n) = 2.0 * yy(m,n)
!
!  Initialize the code.
!
      iflg = 0

      call blktri(iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,
     &  yy,ierror,w)

      if ( ierror .ne. 0 ) then
        write ( *, * ) 'Initial BLKTRI returned IERROR = ',ierror
        ipass = 0
        return
      end if
!
!  Solve the system.
!
      iflg = 1

      call blktri(iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,
     &  yy,ierror,w)

      if ( ierror .ne. 0 ) then
        write ( *, * ) 'Solve BLKTRI returned IERROR = ',ierror
        ipass = 0
        return
      end if
!
!  Compute the discretization error.
!
      errmax = 0.0
      imax = 0
      jmax = 0

      do j = 1, n
        do i = 1, m
          true = x(i)**2 + y(j)**2
          err = abs ( yy(i,j) - true ) 
          if ( err .ge. errmax ) then
            imax = i
            jmax = j
            errmax = err
          end if
        end do
      end do

      write (*,*) ' '
      write (*,*) '  Maximum error is ', errmax
      write (*,*) '  I,J = ',imax, jmax

      return
      end