#include <fintrf.h>
c
c purpose       mexsepeli solves for the 2nd order finite 
c               difference approximation to the separable
c               elliptic equation arising from conformal
c               mapping and grid generation.
c
c
c
c 
c The calling sequence from MATLAB should be
c
c   >> [x, y, perturb, ierror] = mexsepeli (   x, y, ...
c                                              l2, m2, ...
c                                              seta, sxi, ...
c                                              nwrk, nx2 );
c
c
c Output Parameters
c     x,y:
c          represent "x" and "y" grid coordinates
c     perturb:
c          Optional.  Don't think this is needed for our 
c          application.
c     ierror:
c          Optional.  An error flag that indicates invalid 
c          input parameters or failure to find a solution.
c                          = 0 no error
c                          = 4 if attempt to find a solution fails.
c                              (the linear system generated is not
c                              diagonally dominant.)
c
c Input Parameters
c     x,y:
c          represent "x" and "y" grid coordinates
c     l2:
c          x independent variable ranges from 0 to l2
c     m2:
c          y independent variable ranges from 0 to m2
c     seta;
c          digitized points on boundaries 1,3
c     sxi:
c          digitized points on boundaries 2,4
c          
c
c So 
c    1)  nlhs = 2
c    2)  plhs(1) = x
c        plhs(2) = y
c    3)  nrhs = 6 
c    4)  prhs(1) = x 
c        prhs(2) = y 
c        prhs(3) = l2
c        prhs(4) = m2
c        prhs(5) = seta
c        prhs(6) = sxi
c


c
c $Id: mexsepeli.F,v 1.3 2004/09/30 20:57:43 jevans Exp $
c Currently locked by $Locker:  $ (not locked if blank)
c (Time in GMT, EST=GMT-5:00
c $Log: mexsepeli.F,v $
c Revision 1.3  2004/09/30 20:57:43  jevans
c Several automatic casts (from real*8 to complex) were being done that g77
c was complaining about.  Some new complex variables were added that made
c the cast explicit.  g77 no longer spits out any warnings.
c
c Revision 1.2  2004/09/30 20:32:39  jevans
c Removed all declarations of variables that the solaris compiler claimed
c were never used.  Seems to run fine.
c
c Revision 1.1.1.1  2004/09/30 15:22:15  jevans
c initial import
c
c Revision 1.1.1.1  2004/03/02 16:53:17  jevans
c MEXSEPELI
c
c Revision 1.4  2003/02/09 17:51:14  jevans
c Converted some more CMPLX statements to DCMPLX.  Linux compilation didn't
c catch it, but Alpha did.
c
c Revision 1.2  2003/01/15 23:23:04  jevans
c Lots of matrices redefined as for their sizes.  Who knows if it works?
c How the hell did it ever work before?  God how I luv fortran.
c
c Revision 1.1.1.1  2003/01/15 22:31:12  jevans
c Imported sources
c
c
c


c 
c     $Id: mexsepeli.F,v 1.3 2004/09/30 20:57:43 jevans Exp $
c     Currently locked by $Locker:  $ (not locked if blank)
c     (Time in GMT, EST=GMT-5:00
c     $Log: mexsepeli.F,v $
c     Revision 1.3  2004/09/30 20:57:43  jevans
c     Several automatic casts (from real*8 to complex) were being done that g77
c     was complaining about.  Some new complex variables were added that made
c     the cast explicit.  g77 no longer spits out any warnings.
c
c     Revision 1.2  2004/09/30 20:32:39  jevans
c     Removed all declarations of variables that the solaris compiler claimed
c     were never used.  Seems to run fine.
c
c     Revision 1.1.1.1  2004/09/30 15:22:15  jevans
c     initial import
c
c     Revision 1.1.1.1  2004/03/02 16:53:17  jevans
c     MEXSEPELI
c
c     Revision 1.4  2003/02/09 17:51:14  jevans
c     Converted some more CMPLX statements to DCMPLX.  Linux compilation didn't
c     catch it, but Alpha did.
c
c     Revision 1.2  2003/01/15 23:23:04  jevans
c     Lots of matrices redefined as for their sizes.  Who knows if it works?
c     How the hell did it ever work before?  God how I luv fortran.
c
c     Revision 1.1.1.1  2003/01/15 22:31:12  jevans
c     Imported sources
c
c Revision 1.2  1997/10/03  19:43:04  jevans
c Had to change compile time options to include "-align dcommons"
c and "-r8".
c
c the "mexCopyReal8ToPtr" and related functions suffer from the
c fact that they assume all arrays to be a single column.  Since
c my x and y are dimensioned to be larger than their logical
c size (i.e. they are 800x800, but might logically be only 100x80
c or so), those arrays had to be columnified before passing them
c back to MATLAB, and "uncolumnified" upon passing them into
c MATLAB.
c
c Revision 1.1  1997/10/03  17:50:43  jevans
c Initial revision
c
c 
c  




      subroutine mexFunction ( nlhs, plhs, nrhs, prhs )
crps      POINTER plhs(*), prhs(*)
      integer plhs(*), prhs(*)
      integer nlhs, nrhs


      include 'mexsepeli.inc'



c
c     These are the number of rows and columns (m by n) of the
c     matrices we work with.
      integer size_m, size_n, size_mn
      integer i, j, k


c
c     Need some temporary space to store the double-precision
c     matlab arguments.  The correct number of rows and columns
c     is not known at first.
	  real*8 tempx(0:nx2,0:ny2), tempy(0:nx2,0:ny2)


c
c     input ranges for independent variables.
c     0 < x < l2, 0 < y < y2
      integer l2, m2


c
c     Pointers to input arguments
crps      POINTER p 
      integer p 

      real*8 pa(1)
      real*8 PERTRB

c
c     Check for proper number of arguments.
c     There must be at least 1 argument on the left hand size and 
c     no more than 3.
      if ( nlhs .lt. 1 ) then
         call mexErrMsgTxt ( 'Not enough input args for mexsepeli.' )
      elseif ( nlhs .gt. 3 ) then
         call mexErrMsgTxt ( 'Too many input args for mexsepeli.' )
      endif


      if ( nrhs .ne. 6 ) then
         call mexErrMsgTxt('mexsepeli requires 6 input arguments')
      endif


c
c     Check that all input arguments were numeric.
      if ( mxIsNumeric(prhs(1)) .ne. 1 ) then
         call mexErrMsgTxt 
     +        ( 'x input array should be numeric' )
      endif
      if ( mxIsNumeric(prhs(2)) .ne. 1 ) then
         call mexErrMsgTxt 
     +        ( 'y input array should be numeric' )
      endif
      if ( mxIsNumeric(prhs(3)) .ne. 1 ) then
         call mexErrMsgTxt 
     +        ( 'l2 input should be numeric' )
      endif
      if ( mxIsNumeric(prhs(4)) .ne. 1 ) then
         call mexErrMsgTxt 
     +        ( 'm2 input should be numeric' )
      endif
      if ( mxIsNumeric(prhs(5)) .ne. 1 ) then
         call mexErrMsgTxt 
     +        ( 'seta input array should be numeric' )
      endif
      if ( mxIsNumeric(prhs(6)) .ne. 1 ) then
         call mexErrMsgTxt 
     +        ( 'sxi input array should be numeric' )
      endif


c
c     Check to see that the l2 and m2 arguments were scalars
      size_m = mxGetM ( prhs(3) )
      size_n = mxGetN ( prhs(3) )
	  if ( size_n .ne. 1 .or. size_m .ne. 1 ) then
         call mexErrMsgTxt
     +        ('l2 argument should be a scalar' )
      endif

      size_m = mxGetM ( prhs(4) )
      size_n = mxGetN ( prhs(4) )
	  if ( size_n .ne. 1 .or. size_m .ne. 1 ) then
         call mexErrMsgTxt
     +        ('m2 argument should be a scalar' )
      endif




c
c     Load the data into Fortran matrices.

c
c     x
      p = mxGetPr(prhs(1))
      size_m = mxGetM ( prhs(1) )
      size_n = mxGetN ( prhs(1) )
      size_mn = size_m*size_n
      call mxCopyPtrToReal8(p,tempx,size_mn)

      k = 0
      do 10 j = 0, size_n - 1
         do 5 i = 0, size_m - 1
            x(i,j) = tempx(k,0)
            k = k + 1 
    5    continue
   10 continue  
           


c
c     y
      p = mxGetPr(prhs(2))
      size_m = mxGetM ( prhs(2) )
      size_n = mxGetN ( prhs(2) )
      size_mn = size_m*size_n
      call mxCopyPtrToReal8(p,tempy,size_mn)

      k = 0
      do 20 j = 0, size_n - 1
         do 15 i = 0, size_m - 1
            y(i,j) = tempy(k,0)
            k = k + 1 
   15    continue
   20 continue  
           

c
c     l2
      p = mxGetPr(prhs(3))
      call mxCopyPtrToReal8(p,pa,1)
      l2 = pa(1)
      
c
c     m2
      p = mxGetPr(prhs(4))
      call mxCopyPtrToReal8(p,pa,1)
      m2 = pa(1)
      


c
c     seta
      p = mxGetPr(prhs(5))
      size_m = mxGetM ( prhs(5) )
      size_n = mxGetN ( prhs(5) )
      neta = size_m*size_n
      call mxCopyPtrToReal8(p,seta,neta)



c
c     sxi
      p = mxGetPr(prhs(6))
      size_m = mxGetM ( prhs(6) )
      size_n = mxGetN ( prhs(6) )
      nxi = size_m*size_n
      call mxCopyPtrToReal8(p,sxi,nxi)




c
c     Construct right hand side
      do 30 j = 0, neta - 1
        do 25 i = 0, nxi - 1
            rhs(i,j) = 0.
  25    continue
  30  continue


c
c     Call the computational subroutine.
      ewrk(1) = nwrk
      call sepeli(0,2,0.d0,dble(l2),l2,1,wrk,wrk,wrk,wrk,0.d0,
     *    dble(m2),m2,
     *    1,wrk,wrk,wrk,wrk,rhs,x,nx2+1,ewrk,pertrb,ierr)
      if ( ierr .eq. 1 ) then
         call mexErrMsgTxt
     +        ( 'range of independent variables is out of whack' )
      else if ( ierr .eq. 2 ) then
         call mexErrMsgTxt
     +        ( 'boundary condition mbdcnd wrongly specified' )
      else if ( ierr .eq. 3 ) then
         call mexErrMsgTxt
     +        ( 'boundary condition nbdcnd wrongly specified' )
      else if ( ierr .eq. 4 ) then
         call mexErrMsgTxt
     +        ( 'linear system generated is not diagonally dominant' )
      else if ( ierr .eq. 5 ) then
         call mexErrMsgTxt
     +        ( 'idmn is too small' )
      else if ( ierr .eq. 6 ) then
         call mexErrMsgTxt
     +        ( 'm is too small or too large' )
      else if ( ierr .eq. 7 ) then
         call mexErrMsgTxt
     +        ( 'n is too small or too large' )
      else if ( ierr .eq. 8 ) then
         call mexErrMsgTxt
     +        ( 'iorder is not 2 or 4' )
      else if ( ierr .eq. 9 ) then
         call mexErrMsgTxt
     +        ( 'intl is not 0 or 1' )
      else if ( ierr .eq. 10 ) then
         call mexErrMsgTxt
     +        ( 'afun*dfun less than or equal to 0' )
      else if ( ierr .eq. 11 ) then
         call mexErrMsgTxt
     +        ( 'work space length input in w(1) is not right' )
      end if
          



      ewrk(1) = nwrk
      call sepeli(0,2,0.d0,dble(l2),l2,1,wrk,wrk,wrk,wrk,0.d0,
     *    dble(m2),m2,
     *    1,wrk,wrk,wrk,wrk,rhs,y,nx2+1,ewrk,pertrb,ierr)
      if ( ierr .eq. 1 ) then
         call mexErrMsgTxt
     +        ( 'range of independent variables is out of whack' )
      else if ( ierr .eq. 2 ) then
         call mexErrMsgTxt
     +        ( 'boundary condition mbdcnd wrongly specified' )
      else if ( ierr .eq. 3 ) then
         call mexErrMsgTxt
     +        ( 'boundary condition nbdcnd wrongly specified' )
      else if ( ierr .eq. 4 ) then
         call mexErrMsgTxt
     +        ( 'linear system generated is not diagonally dominant' )
      else if ( ierr .eq. 5 ) then
         call mexErrMsgTxt
     +        ( 'idmn is too small' )
      else if ( ierr .eq. 6 ) then
         call mexErrMsgTxt
     +        ( 'm is too small or too large' )
      else if ( ierr .eq. 7 ) then
         call mexErrMsgTxt
     +        ( 'n is too small or too large' )
      else if ( ierr .eq. 8 ) then
         call mexErrMsgTxt
     +        ( 'iorder is not 2 or 4' )
      else if ( ierr .eq. 9 ) then
         call mexErrMsgTxt
     +        ( 'intl is not 0 or 1' )
      else if ( ierr .eq. 10 ) then
         call mexErrMsgTxt
     +        ( 'afun*dfun less than or equal to 0' )
      else if ( ierr .eq. 11 ) then
         call mexErrMsgTxt
     +        ( 'work space length input in w(1) is not right' )
      end if




c
c     Create matrices for output arguments
c     Since the mxCopyReal8ToPtr routine just assumes one long
c     column matrix, we have to "columnify" both x and y.

c
c     x
      size_m = mxGetM ( prhs(1) )
      size_n = mxGetN ( prhs(1) )
      size_mn = size_m*size_n
      plhs(1) = mxCreateDoubleMatrix(size_m, size_n, 0)
      p = mxGetPr(plhs(1))

      k = 0
      do 40 j = 0, size_n - 1
         do 35 i = 0, size_m - 1
            tempx(k,0) = x(i,j)
            k = k + 1 
   35    continue
   40 continue  
           
      call mxCopyReal8ToPtr ( tempx, p, size_mn )


c
c     y
      size_m = mxGetM ( prhs(2) )
      size_n = mxGetN ( prhs(2) )
      size_mn = size_m*size_n
      plhs(2) = mxCreateDoubleMatrix(size_m, size_n, 0)
      p = mxGetPr(plhs(2))

      k = 0
      do 50 j = 0, size_n - 1
         do 45 i = 0, size_m - 1
            tempy(k,0) = y(i,j)
            k = k + 1 
   45    continue
   50 continue  
           
      call mxCopyReal8ToPtr ( tempy, p, size_mn )



      return
      end





c      Subroutine SEPELI(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,
c     *    N,NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,GRHS,USOL,IDMN,W,PERTRB,
c     *    IERROR)
      Subroutine SEPELI(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,
     *    N,NBDCND,BDC,GAMA,BDD,XNU,GRHS,USOL,IDMN,W,PERTRB,
     *    IERROR)
c
c dimension of           bda(n+1), bdb(n+1), bdc(m+1), bdd(m+1),
c arguments              usol(idmn,n+1),grhs(idmn,n+1),
c                        w (see argument list)
c
c latest revision        march 1985
c
c purpose                sepeli solves for either the second-order
c                        finite difference approximation or a
c                        fourth-order approximation to a separable
c                        elliptic equation
c
c                          2    2
c                          af(x)*d u/dx + bf(x)*du/dx  + cf(x)*u +
c                          2    2
c                          df(y)*d u/dy  + ef(y)*du/dy + ff(y)*u
c
c                          = g(x,y)
c
c                        on a rectangle (x greater than or equal to a
c                        and less than or equal to b; y greater than
c                        or equal to c and less than or equal to d).
c                        any combination of periodic or mixed boundary
c                        conditions is allowed.
c
c                        the possible boundary conditions are:
c                        in the x-direction:
c                        (0) periodic, u(x+b-a,y)=u(x,y) for all
c                            y,x (1) u(a,y), u(b,y) are specified for
c                            all y
c                        (2) u(a,y), du(b,y)/dx+beta*u(b,y) are
c                            specified for all y
c                        (3) du(a,y)/dx+alpha*u(a,y),du(b,y)/dx+
c                            beta*u(b,y) are specified for all y
c                        (4) du(a,y)/dx+alpha*u(a,y),u(b,y) are
c                            specified for all y
c
c                        in the y-direction:
c                        (0) periodic, u(x,y+d-c)=u(x,y) for all x,y
c                        (1) u(x,c),u(x,d) are specified for all x
c                        (2) u(x,c),du(x,d)/dy+xnu*u(x,d) are
c                            specified for all x
c                        (3) du(x,c)/dy+gama*u(x,c),du(x,d)/dy+
c                            xnu*u(x,d) are specified for all x
c                        (4) du(x,c)/dy+gama*u(x,c),u(x,d) are
c                            specified for all x
c
c usage                  call sepeli (intl,iorder,a,b,m,mbdcnd,bda,
c                                     alpha,bdb,beta,c,d,n,nbdcnd,bdc,
c                                     gama,bdd,xnu,cofx,cofy,grhs,usol,
c                                     idmn,w,pertrb,ierror)
c
c arguments
c on input               intl
c                          = 0 on initial entry to sepeli or if any
c                              of the arguments c,d, n, nbdcnd, cofy
c                              are changed from a previous call
c                          = 1 if c, d, n, nbdcnd, cofy are unchanged
c                              from the previous call.
c
c                        iorder
c                          = 2 if a second-order approximation
c                              is sought
c                          = 4 if a fourth-order approximation
c                              is sought
c
c                        a,b
c                          the range of the x-independent variable,
c                          i.e., x is greater than or equal to a
c                          and less than or equal to b.  a must be
c                          less than b.
c
c                        m
c                          the number of panels into which the
c                          interval [a,b] is subdivided. hence,
c                          there will be m+1 grid points in the x-
c                          direction given by xi=a+(i-1)*dlx
c                          for i=1,2,...,m+1 where dlx=(b-a)/m is
c                          the panel width.  m must be less than
c                          idmn and greater than 5.
c
c                        mbdcnd
c                          indicates the type of boundary condition
c                          at x=a and x=b
c
c                          = 0 if the solution is periodic in x, i.e.,
c                              u(x+b-a,y)=u(x,y)  for all y,x
c                          = 1 if the solution is specified at x=a
c                              and x=b, i.e., u(a,y) and u(b,y) are
c                              specified for all y
c                          = 2 if the solution is specified at x=a and
c                              the boundary condition is mixed at x=b,
c                              i.e., u(a,y) and du(b,y)/dx+beta*u(b,y)
c                              are specified for all y
c                          = 3 if the boundary conditions at x=a and
c                              x=b are mixed, i.e.,
c                              du(a,y)/dx+alpha*u(a,y) and
c                              du(b,y)/dx+beta*u(b,y) are specified
c                              for all y
c                          = 4 if the boundary condition at x=a is
c                              mixed and the solution is specified
c                              at x=b, i.e., du(a,y)/dx+alpha*u(a,y)
c                              and u(b,y) are specified for all y
c
c                        bda
c                          a one-dimensional array of length n+1
c                          that specifies the values of
c                          du(a,y)/dx+ alpha*u(a,y) at x=a, when
c                          mbdcnd=3 or 4.
c                          bda(j) = du(a,yj)/dx+alpha*u(a,yj),
c                          j=1,2,...,n+1. when mbdcnd has any other
c                          other value, bda is a dummy parameter.
c
c                        alpha
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at x=a
c                          (see argument bda).  if mbdcnd is not
c                          equal to 3 or 4 then alpha is a dummy
c                          parameter.
c
c                        bdb
c                          a one-dimensional array of length n+1
c                          that specifies the values of
c                          du(b,y)/dx+ beta*u(b,y) at x=b.
c                          when mbdcnd=2 or 3
c                          bdb(j) = du(b,yj)/dx+beta*u(b,yj),
c                          j=1,2,...,n+1. when mbdcnd has any other
c                          other value, bdb is a dummy parameter.
c
c                        beta
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at
c                          x=b (see argument bdb).  if mbdcnd is
c                          not equal to 2 or 3 then beta is a dummy
c                          parameter.
c
c                        c,d
c                          the range of the y-independent variable,
c                          i.e., y is greater than or equal to c
c                          and less than or equal to d.  c must be
c                          less than d.
c
c                        n
c                          the number of panels into which the
c                          interval [c,d] is subdivided.
c                          hence, there will be n+1 grid points
c                          in the y-direction given by
c                          yj=c+(j-1)*dly for j=1,2,...,n+1 where
c                          dly=(d-c)/n is the panel width.
c                          in addition, n must be greater than 4.
c
c                        nbdcnd
c                          indicates the types of boundary conditions
c                          at y=c and y=d
c
c                          = 0 if the solution is periodic in y,
c                              i.e., u(x,y+d-c)=u(x,y)  for all x,y
c                          = 1 if the solution is specified at y=c
c                              and y = d, i.e., u(x,c) and u(x,d)
c                              are specified for all x
c                          = 2 if the solution is specified at y=c
c                              and the boundary condition is mixed
c                              at y=d, i.e., u(x,c) and
c                              du(x,d)/dy+xnu*u(x,d) are specified
c                              for all x
c                          = 3 if the boundary conditions are mixed
c                              at y=c and y=d, i.e.,
c                              du(x,d)/dy+gama*u(x,c) and
c                              du(x,d)/dy+xnu*u(x,d) are specified
c                              for all x
c                          = 4 if the boundary condition is mixed
c                              at y=c and the solution is specified
c                              at y=d, i.e. du(x,c)/dy+gama*u(x,c)
c                              and u(x,d) are specified for all x
c
c                        bdc
c                          a one-dimensional array of length m+1
c                          that specifies the value of
c                          du(x,c)/dy+gama*u(x,c) at y=c.
c                          when nbdcnd=3 or 4 bdc(i) = du(xi,c)/dy +
c                          gama*u(xi,c), i=1,2,...,m+1.
c                          when nbdcnd has any other value, bdc
c                          is a dummy parameter.
c
c                        gama
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at
c                          y=c (see argument bdc).  if nbdcnd is
c                          not equal to 3 or 4 then gama is a dummy
c                          parameter.
c
c                        bdd
c                          a one-dimensional array of length m+1
c                          that specifies the value of
c                          du(x,d)/dy + xnu*u(x,d) at y=c.
c                          when nbdcnd=2 or 3 bdd(i) = du(xi,d)/dy +
c                          xnu*u(xi,d), i=1,2,...,m+1.
c                          when nbdcnd has any other value, bdd
c                          is a dummy parameter.
c
c                        xnu
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at
c                          y=d (see argument bdd).  if nbdcnd is
c                          not equal to 2 or 3 then xnu is a
c                          dummy parameter.
c
c                        cofx
c                          a user-supplied subprogram with
c                          parameters x, afun, bfun, cfun which
c                          returns the values of the x-dependent
c                          coefficients af(x), bf(x), cf(x) in the
c                          elliptic equation at x.
c
c                        cofy
c                          a user-supplied subprogram with parameters
c                          y, dfun, efun, ffun which returns the
c                          values of the y-dependent coefficients
c                          df(y), ef(y), ff(y) in the elliptic
c                          equation at y.
c
c                          note:  cofx and cofy must be declared
c                          external in the calling routine.
c                          the values returned in afun and dfun
c                          must satisfy afun*dfun greater than 0
c                          for a less than x less than b, c less
c                          than y less than d (see ierror=10).
c                          the coefficients provided may lead to a
c                          matrix equation which is not diagonally
c                          dominant in which case solution may fail
c                          (see ierror=4).
c
c                        grhs
c                          a two-dimensional array that specifies the
c                          values of the right-hand side of the
c                          elliptic equation, i.e.,
c                          grhs(i,j)=g(xi,yi), for i=2,...,m,
c                          j=2,...,n.  at the boundaries, grhs is
c                          defined by
c
c                          mbdcnd   grhs(1,j)   grhs(m+1,j)
c                          ------   ---------   -----------
c                            0      g(a,yj)     g(b,yj)
c                            1         *           *
c                            2         *        g(b,yj)  j=1,2,...,n+1
c                            3      g(a,yj)     g(b,yj)
c                            4      g(a,yj)        *
c
c                          nbdcnd   grhs(i,1)   grhs(i,n+1)
c                          ------   ---------   -----------
c                            0      g(xi,c)     g(xi,d)
c                            1         *           *
c                            2         *        g(xi,d)  i=1,2,...,m+1
c                            3      g(xi,c)     g(xi,d)
c                            4      g(xi,c)        *
c
c                          where * means these quantities are not used.
c                          grhs should be dimensioned idmn by at least
c                          n+1 in the calling routine.
c
c                        usol
c                          a two-dimensional array that specifies the
c                          values of the solution along the boundaries.
c                          at the boundaries, usol is defined by
c
c                          mbdcnd   usol(1,j)   usol(m+1,j)
c                          ------   ---------   -----------
c                            0         *           *
c                            1      u(a,yj)     u(b,yj)
c                            2      u(a,yj)        *     j=1,2,...,n+1
c                            3         *           *
c                            4         *        u(b,yj)
c
c                          nbdcnd   usol(i,1)   usol(i,n+1)
c                          ------   ---------   -----------
c                            0         *           *
c                            1      u(xi,c)     u(xi,d)
c                            2      u(xi,c)        *     i=1,2,...,m+1
c                            3         *           *
c                            4         *        u(xi,d)
c
c                          where * means the quantities are not used
c                          in the solution.
c
c                          if iorder=2, the user may equivalence grhs
c                          and usol to save space.  note that in this
c                          case the tables specifying the boundaries
c                          of the grhs and usol arrays determine the
c                          boundaries uniquely except at the corners.
c                          if the tables call for both g(x,y) and
c                          u(x,y) at a corner then the solution must
c                          be chosen.  for example, if mbdcnd=2 and
c                          nbdcnd=4, then u(a,c), u(a,d), u(b,d) must
c                          be chosen at the corners in addition
c                          to g(b,c).
c
c                          if iorder=4, then the two arrays, usol and
c                          grhs, must be distinct.
c
c                          usol should be dimensioned idmn by at least
c                          n+1 in the calling routine.
c
c                        idmn
c                          the row (or first) dimension of the arrays
c                          grhs and usol as it appears in the program
c                          calling sepeli.  this parameter is used
c                          to specify the variable dimension of grhs
c                          and usol.  idmn must be at least 7 and
c                          greater than or equal to m+1.
c
c                        w
c                          a one-dimensional array that must be
c                          provided by the user for work space.
c                          let k=int(log2(n+1))+1 and set l=2**(k+1).
c                          then (k-2)*l+k+10*n+12*m+27 will suffice
c                          as a length of w.  the actual length of w
c                          in the calling routine must be set in w(1)
c                          (see ierror=11).
c
c on output              usol
c                          contains the approximate solution to the
c                          elliptic equation.
c                          usol(i,j) is the approximation to u(xi,yj)
c                          for i=1,2...,m+1 and j=1,2,...,n+1.
c                          the approximation has error
c                          o(dlx**2+dly**2) if called with iorder=2
c                          and o(dlx**4+dly**4) if called with
c                          iorder=4.
c
c                        w
c                          contains intermediate values that must not
c                          be destroyed if sepeli is called again with
c                          intl=1.  in addition w(1) contains the
c                          exact minimal length (in floating point)
c                          required for the work space (see ierror=11).
c
c                        pertrb
c                          if a combination of periodic or derivative
c                          boundary conditions
c                          (i.e., alpha=beta=0 if mbdcnd=3;
c                          gama=xnu=0 if nbdcnd=3) is specified
c                          and if the coefficients of u(x,y) in the
c                          separable elliptic equation are zero
c                          (i.e., cf(x)=0 for x greater than or equal
c                          to a and less than or equal to b;
c                          ff(y)=0 for y greater than or equal to c
c                          and less than or equal to d) then a
c                          solution may not exist.  pertrb is a
c                          constant calculated and subtracted from
c                          the right-hand side of the matrix equations
c                          generated by sepeli which insures that a
c                          solution exists. sepeli then computes this
c                          solution which is a weighted minimal least
c                          squares solution to the original problem.
c
c                        ierror
c                          an error flag that indicates invalid input
c                          parameters or failure to find a solution
c                          = 0 no error
c                          = 1 if a greater than b or c greater than d
c                          = 2 if mbdcnd less than 0 or mbdcnd greater
c                              than 4
c                          = 3 if nbdcnd less than 0 or nbdcnd greater
c                              than 4
c                          = 4 if attempt to find a solution fails.
c                              (the linear system generated is not
c                              diagonally dominant.)
c                          = 5 if idmn is too small
c                              (see discussion of idmn)
c                          = 6 if m is too small or too large
c                              (see discussion of m)
c                          = 7 if n is too small (see discussion of n)
c                          = 8 if iorder is not 2 or 4
c                          = 9 if intl is not 0 or 1
c                          = 10 if afun*dfun less than or equal to 0
c                               for some interior mesh point (xi,yj)
c                          = 11 if the work space length input in w(1)
c                               is less than the exact minimal work
c                               space length required output in w(1).
c
c                          note (concerning ierror=4):  for the
c                          coefficients input through cofx, cofy,
c                          the discretization may lead to a block
c                          tridiagonal linear system which is not
c                          diagonally dominant (for example, this
c                          happens if cfun=0 and bfun/(2.*dlx) greater
c                          than afun/dlx**2).  in this case solution
c                          may fail.  this cannot happen in the limit
c                          as dlx, dly approach zero.  hence, the
c                          condition may be remedied by taking larger
c                          values for m or n.
c
c special conditions     see cofx, cofy argument descriptions above.
c
c i/o                    none
c
c precision              single
c
c required library       blktri, comf, q8qst4, and sepaux, which are
c files                  loaded by default on ncar's cray machines.
c
c language               fortran
c
c history                developed at ncar during 1975-76 by
c                        john c. adams of the scientific computing
c                        division.  released on ncar's public software
c                        libraries in january 1980.
c
c portability            fortran 66
c
c algorithm              sepeli automatically discretizes the
c                        separable elliptic equation which is then
c                        solved by a generalized cyclic reduction
c                        algorithm in the subroutine, blktri.  the
c                        fourth-order solution is obtained using
c                        'deferred corrections' which is described
c                        and referenced in sections, references and
c                        method.
c
c timing                 the operational count is proportional to
c                        m*n*log2(n).
c
c accuracy               the following accuracy results were obtained
c                        on a cdc 7600.  note that the fourth-order
c                        accuracy is not realized until the mesh is
c                        sufficiently refined.
c
c                                     second-order  fourth-order
c                            m    n     error         error
c
c                             6    6    6.8e-1        1.2e0
c                            14   14    1.4e-1        1.8e-1
c                            30   30    3.2e-2        9.7e-3
c                            62   62    7.5e-3        3.0e-4
c                           126  126    1.8e-3        3.5e-6
c
c portability            fortran 66
c
c references             keller, h.b., numerical methods for two-point
c                        boundary-value problems, blaisdel (1968),
c                        waltham, mass.
c
c                        swarztrauber, p., and r. sweet (1975):
c                        efficient fortran subprograms for the
c                        solution of elliptic partial differential
c                        equations.  ncar technical note
c                        ncar-tn/ia-109, pp. 135-137.
c***********************************************************************
      IMPLICIT NONE
      integer IDMN, INTL, IORDER, MBDCND, NBDCND, M, N, IERROR
      integer L, LOGB2N, LL, K, LENGTH, LINPUT
      integer LOUTPT
      integer I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11, I12, I13
      Dimension GRHS(IDMN,1), USOL(IDMN,IDMN)
      Dimension BDA(1), BDB(1), BDC(1), BDD(1), W(1)
      real*8 BDA
      real*8 ALPHA
      real*8 BDB
      real*8 BETA
      real*8 BDC
      real*8 GAMA
      real*8 BDD
      real*8 XNU
      real*8 GRHS
      real*8 USOL
      real*8 W
      real*8 A, B, C, D, PERTRB
c      External COFX, COFY
      Logical Q8Q4
      Save Q8Q4
      Data Q8Q4 /.TRUE./
c
      If (Q8Q4) Then
c         call q8qst4('loclib','sepeli','sepeli','version 01')
        Q8Q4 = .FALSE.
      End If
c
c     check input parameters
c
c      Call CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,COFX,COFY,IDMN,
c     *    IERROR)
      Call CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,IDMN,
     *    IERROR)
      If (IERROR.NE.0) Return
c
c     compute minimum work space and check work space length input
c
      L = N + 1
      If (NBDCND.EQ.0) L = N
      LOGB2N = INT(DLOG(DBLE(L)+0.d5)/DLOG(2.d0)) + 1
      LL = 2 ** (LOGB2N+1)
      K = M + 1
      L = N + 1
      LENGTH = (LOGB2N-2) * LL + LOGB2N + MAX0(2*L,6*K) + 5
      If (NBDCND.EQ.0) LENGTH = LENGTH + 2 * L
      IERROR = 11
      LINPUT = INT(W(1)+0.5)
      LOUTPT = LENGTH + 6 * (K+L) + 1
      W(1) = DBLE(LOUTPT)
      If (LOUTPT.GT.LINPUT) Return
      IERROR = 0
c
c     set work space indices
c
      I1 = LENGTH + 2
      I2 = I1 + L
      I3 = I2 + L
      I4 = I3 + L
      I5 = I4 + L
      I6 = I5 + L
      I7 = I6 + L
      I8 = I7 + K
      I9 = I8 + K
      I10 = I9 + K
      I11 = I10 + K
      I12 = I11 + K
      I13 = 2
c      Call SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N,
c     *    NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,W(I1),W(I2),W(I3),W(I4),
c     *    W(I5),W(I6),W(I7),W(I8),W(I9),W(I10),W(I11),W(I12),GRHS,USOL,
c     *    IDMN,W(I13),PERTRB,IERROR)
      Call SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N,
     *    NBDCND,BDC,GAMA,BDD,XNU,W(I1),W(I2),W(I3),W(I4),
     *    W(I5),W(I6),W(I7),W(I8),W(I9),W(I10),W(I11),W(I12),GRHS,USOL,
     *    IDMN,W(I13),PERTRB,IERROR)
      Return
      End













c      Subroutine SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,
c     *    N,NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,AN,BN,CN,DN,UN,ZN,AM,BM,
c     *    CM,DM,UM,ZM,GRHS,USOL,IDMN,W,PERTRB,IERROR)
      Subroutine SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,
     *    N,NBDCND,BDC,GAMA,BDD,XNU,AN,BN,CN,DN,UN,ZN,AM,BM,
     *    CM,DM,UM,ZM,GRHS,USOL,IDMN,W,PERTRB,IERROR)
c
c     spelip sets up vectors and arrays for input to blktri
c     and computes a second order solution in usol.  a return jump to
c     sepeli occurrs if iorder=2.  if iorder=4 a fourth order
c     solution is generated in usol.
c
      IMPLICIT NONE
      integer IDMN
      Dimension BDA(1), BDB(1), BDC(1), BDD(1), W(1)
      Dimension GRHS(IDMN,1), USOL(IDMN,IDMN)
      Dimension AN(1), BN(1), CN(1), DN(1), UN(1), ZN(1)
      Dimension AM(1), BM(1), CM(1), DM(1), UM(1), ZM(1)

      integer KSWX, KSWY, K, MBDCND, NBDCND 

      integer AIT, BIT, CIT, DIT, IS, MIT, NIT, L
      real*8 MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4
      Common /SPLP/ KSWX, KSWY, K, L, 
     *    MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4,
     *    AIT, BIT, CIT, DIT, MIT, NIT, IS

      Logical SINGLR

      real*8 BDA, BDB, BDC, BDD
      real*8 ALPHA, BETA, GAMA
      real*8 A, B, C, D, W
      real*8 XNU
      real*8 GRHS, USOL
      real*8 AN, BN, CN, DN, UN, ZN
      real*8 AM, BM, CM, DM, UM, ZM
      real*8 YJ, DJ, EJ, FJ
      real*8 XI, AI, BI, CI
      real*8 PRTRB, PERTRB

      integer M, N, I, J, I1, MP, NP, IORDER, IORD, INTL, IERROR
      real*8 AXI, BXI, CXI, DYJ, EYJ, FYJ, AX1, CXM, DY1, FYN

c      External COFX, COFY
c
c     set parameters internally
c
      KSWX = MBDCND + 1
      KSWY = NBDCND + 1
      K = M + 1
      L = N + 1
      AIT = A
      BIT = B
      CIT = C
      DIT = D
c
c     set right hand side values from grhs in usol on the interior
c     and non-specified boundaries.
c
      Do 20 I = 2, M
        Do 10 J = 2, N
          USOL(I,J) = GRHS(I,J)
   10   Continue
   20 Continue
      If (KSWX.NE.2.AND.KSWX.NE.3) Then
        Do 30 J = 2, N
          USOL(1,J) = GRHS(1,J)
   30   Continue
      End If
      If (KSWX.NE.2.AND.KSWX.NE.5) Then
        Do 40 J = 2, N
          USOL(K,J) = GRHS(K,J)
   40   Continue
      End If
      If (KSWY.NE.2.AND.KSWY.NE.3) Then
        Do 50 I = 2, M
          USOL(I,1) = GRHS(I,1)
   50   Continue
      End If
      If (KSWY.NE.2.AND.KSWY.NE.5) Then
        Do 60 I = 2, M
          USOL(I,L) = GRHS(I,L)
   60   Continue
      End If
      If (KSWX.NE.2.AND.KSWX.NE.3.AND.KSWY.NE.2.AND.KSWY.NE.3) 
     *    USOL(1,1) = GRHS(1,1)
      If (KSWX.NE.2.AND.KSWX.NE.5.AND.KSWY.NE.2.AND.KSWY.NE.3) 
     *    USOL(K,1) = GRHS(K,1)
      If (KSWX.NE.2.AND.KSWX.NE.3.AND.KSWY.NE.2.AND.KSWY.NE.5) 
     *    USOL(1,L) = GRHS(1,L)
      If (KSWX.NE.2.AND.KSWX.NE.5.AND.KSWY.NE.2.AND.KSWY.NE.5) 
     *    USOL(K,L) = GRHS(K,L)
      I1 = 1
c
c     set switches for periodic or non-periodic boundaries
c
      MP = 1
      NP = 1
      If (KSWX.EQ.1) MP = 0
      If (KSWY.EQ.1) NP = 0
c
c     set dlx,dly and size of block tri-diagonal system generated
c     in nint,mint
c
      DLX = (BIT-AIT) / DBLE(M)
      MIT = K - 1
      If (KSWX.EQ.2) MIT = K - 2
      If (KSWX.EQ.4) MIT = K
      DLY = (DIT-CIT) / DBLE(N)
      NIT = L - 1
      If (KSWY.EQ.2) NIT = L - 2
      If (KSWY.EQ.4) NIT = L
      TDLX3 = 2.0 * DLX ** 3
      DLX4 = DLX ** 4
      TDLY3 = 2.0 * DLY ** 3
      DLY4 = DLY ** 4
c
c     set subscript limits for portion of array to input to blktri
c
      IS = 1
      JS = 1
      If (KSWX.EQ.2.OR.KSWX.EQ.3) IS = 2
      If (KSWY.EQ.2.OR.KSWY.EQ.3) JS = 2
      NS = NIT + JS - 1
      MS = MIT + IS - 1
c
c     set x - direction
c
      Do 70 I = 1, MIT
        XI = AIT + DBLE(IS+I-2) * DLX
        Call COFX(XI,AI,BI,CI)
        AXI = (AI/DLX-0.5*BI) / DLX
        BXI = -2. * AI / DLX ** 2 + CI
        CXI = (AI/DLX+0.5*BI) / DLX
        AM(I) = AXI
        BM(I) = BXI
        CM(I) = CXI
   70 Continue
c
c     set y direction
c
      Do 80 J = 1, NIT
        YJ = CIT + DBLE(JS+J-2) * DLY
        Call COFY(YJ,DJ,EJ,FJ)
        DYJ = (DJ/DLY-0.5*EJ) / DLY
        EYJ = (-2.*DJ/DLY**2+FJ)
        FYJ = (DJ/DLY+0.5*EJ) / DLY
        AN(J) = DYJ
        BN(J) = EYJ
        CN(J) = FYJ
   80 Continue
c
c     adjust edges in x direction unless periodic
c
      AX1 = AM(1)
      CXM = CM(MIT)
      Go To (130,90,110,120,100), KSWX
c
c     dirichlet-dirichlet in x direction
c
   90 AM(1) = 0.0
      CM(MIT) = 0.0
      Go To 130
c
c     mixed-dirichlet in x direction
c
  100 AM(1) = 0.0
      BM(1) = BM(1) + 2. * ALPHA * DLX * AX1
      CM(1) = CM(1) + AX1
      CM(MIT) = 0.0
      Go To 130
c
c     dirichlet-mixed in x direction
c
  110 AM(1) = 0.0
      AM(MIT) = AM(MIT) + CXM
      BM(MIT) = BM(MIT) - 2. * BETA * DLX * CXM
      CM(MIT) = 0.0
      Go To 130
c
c     mixed - mixed in x direction
c
  120 Continue
      AM(1) = 0.0
      BM(1) = BM(1) + 2. * DLX * ALPHA * AX1
      CM(1) = CM(1) + AX1
      AM(MIT) = AM(MIT) + CXM
      BM(MIT) = BM(MIT) - 2. * DLX * BETA * CXM
      CM(MIT) = 0.0
  130 Continue
c
c     adjust in y direction unless periodic
c
      DY1 = AN(1)
      FYN = CN(NIT)
      Go To (180,140,160,170,150), KSWY
c
c     dirichlet-dirichlet in y direction
c
  140 Continue
      AN(1) = 0.0
      CN(NIT) = 0.0
      Go To 180
c
c     mixed-dirichlet in y direction
c
  150 Continue
      AN(1) = 0.0
      BN(1) = BN(1) + 2. * DLY * GAMA * DY1
      CN(1) = CN(1) + DY1
      CN(NIT) = 0.0
      Go To 180
c
c     dirichlet-mixed in y direction
c
  160 AN(1) = 0.0
      AN(NIT) = AN(NIT) + FYN
      BN(NIT) = BN(NIT) - 2. * DLY * XNU * FYN
      CN(NIT) = 0.0
      Go To 180
c
c     mixed - mixed direction in y direction
c
  170 Continue
      AN(1) = 0.0
      BN(1) = BN(1) + 2. * DLY * GAMA * DY1
      CN(1) = CN(1) + DY1
      AN(NIT) = AN(NIT) + FYN
      BN(NIT) = BN(NIT) - 2.0 * DLY * XNU * FYN
      CN(NIT) = 0.0
  180 If (KSWX.NE.1) Then
c
c     adjust usol along x edge
c
        Do 190 J = JS, NS
          If (KSWX.EQ.2.OR.KSWX.EQ.3) Then
            USOL(IS,J) = USOL(IS,J) - AX1 * USOL(1,J)
          Else
            USOL(IS,J) = USOL(IS,J) + 2.0 * DLX * AX1 * BDA(J)
          End If
          If (KSWX.EQ.2.OR.KSWX.EQ.5) Then
            USOL(MS,J) = USOL(MS,J) - CXM * USOL(K,J)
          Else
            USOL(MS,J) = USOL(MS,J) - 2.0 * DLX * CXM * BDB(J)
          End If
  190   Continue
      End If
      If (KSWY.NE.1) Then
c
c     adjust usol along y edge
c
        Do 200 I = IS, MS
          If (KSWY.EQ.2.OR.KSWY.EQ.3) Then
            USOL(I,JS) = USOL(I,JS) - DY1 * USOL(I,1)
          Else
            USOL(I,JS) = USOL(I,JS) + 2.0 * DLY * DY1 * BDC(I)
          End If
          If (KSWY.EQ.2.OR.KSWY.EQ.5) Then
            USOL(I,NS) = USOL(I,NS) - FYN * USOL(I,L)
          Else
            USOL(I,NS) = USOL(I,NS) - 2.0 * DLY * FYN * BDD(I)
          End If
  200   Continue
      End If
c
c     save adjusted edges in grhs if iorder=4
c
      If (IORDER.EQ.4) Then
        Do 210 J = JS, NS
          GRHS(IS,J) = USOL(IS,J)
          GRHS(MS,J) = USOL(MS,J)
  210   Continue
        Do 220 I = IS, MS
          GRHS(I,JS) = USOL(I,JS)
          GRHS(I,NS) = USOL(I,NS)
  220   Continue
      End If
      IORD = IORDER
      PERTRB = 0.0
c
c     check if operator is singular
c
c      Call CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,COFX,COFY,SINGLR)
      Call CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,SINGLR)
c
c     compute non-zero eigenvector in null space of transpose
c     if singular
c
      If (SINGLR) Call SEPTRI(MIT,AM,BM,CM,DM,UM,ZM)
      If (SINGLR) Call SEPTRI(NIT,AN,BN,CN,DN,UN,ZN)
c
c     make initialization call to blktri
c
      If (INTL.EQ.0) Call BLKTRI(INTL,NP,NIT,AN,BN,CN,MP,MIT,AM,BM,CM,
     *    IDMN,USOL(IS,JS),IERROR,W)
      If (IERROR.NE.0) Return
c
c     adjust right hand side if necessary
c
  230 Continue
      If (SINGLR) Call SEPORT(USOL,IDMN,ZN,ZM,PERTRB)
c
c     compute solution
c
      Call BLKTRI(I1,NP,NIT,AN,BN,CN,MP,MIT,AM,BM,CM,IDMN,USOL(IS,JS),
     *    IERROR,W)
      If (IERROR.NE.0) Return
c
c     set periodic boundaries if necessary
c
      If (KSWX.EQ.1) Then
        Do 240 J = 1, L
          USOL(K,J) = USOL(1,J)
  240   Continue
      End If
      If (KSWY.EQ.1) Then
        Do 250 I = 1, K
          USOL(I,L) = USOL(I,1)
  250   Continue
      End If
c
c     minimize solution with respect to weighted least squares
c     norm if operator is singular
c
      If (SINGLR) Call SEPMIN(USOL,IDMN,ZN,ZM,PRTRB)
c
c     return if deferred corrections and a fourth order solution are
c     not flagged
c
      If (IORD.EQ.2) Return
      IORD = 2
c
c     compute new right hand side for fourth order solution
c
c      Call DEFER(COFX,COFY,IDMN,USOL,GRHS)
      Call DEFER(IDMN,USOL,GRHS)
      Go To 230
      End



















c      Subroutine CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,COFX,COFY,
c     *    IDMN,IERROR)
      Subroutine CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,
     *    IDMN,IERROR)
c
c     this program checks the input parameters for errors
c
c      External COFX, COFY
c
c     check definition of solution region
c
      IMPLICIT NONE
      integer IERROR, MBDCND, NBDCND, IDMN, M, N, IORDER, INTL
      real*8 A, B, C, D

      real*8 YJ, DJ, EJ, FJ
      real*8 XI, AI, BI, CI, DLX, DLY
      integer I, J




      IERROR = 1
      If (A.GE.B.OR.C.GE.D) Return
c
c     check boundary switches
c
      IERROR = 2
      If (MBDCND.LT.0.OR.MBDCND.GT.4) Return
      IERROR = 3
      If (NBDCND.LT.0.OR.NBDCND.GT.4) Return
c
c     check first dimension in calling routine
c
      IERROR = 5
      If (IDMN.LT.7) Return
c
c     check m
c
      IERROR = 6
      If (M.GT.(IDMN-1).OR.M.LT.6) Return
c
c     check n
c
      IERROR = 7
      If (N.LT.5) Return
c
c     check iorder
c
      IERROR = 8
      If (IORDER.NE.2.AND.IORDER.NE.4) Return
c
c     check intl
c
      IERROR = 9
      If (INTL.NE.0.AND.INTL.NE.1) Return
c
c     check that equation is elliptic
c
      DLX = (B-A) / DBLE(M)
      DLY = (D-C) / DBLE(N)
      Do 20 I = 2, M
        XI = A + DBLE(I-1) * DLX
        Call COFX(XI,AI,BI,CI)
        Do 10 J = 2, N
          YJ = C + DBLE(J-1) * DLY
          Call COFY(YJ,DJ,EJ,FJ)
          If (AI*DJ.LE.0.0) Then
            IERROR = 10
            Return
          End If
   10   Continue
   20 Continue
c
c     no error found
c
      IERROR = 0
      Return
      End















c      Subroutine CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,COFX,COFY,
c     *    SINGLR)
      Subroutine CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,
     *    SINGLR)
c
c     this subroutine checks if the pde   sepeli
c     must solve is a singular operator
c
      IMPLICIT NONE

      integer MBDCND, NBDCND

      integer KSWX, KSWY, K, L
      integer AIT, BIT, CIT, DIT, IS, MIT, NIT
      real*8 MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4
      Common /SPLP/ KSWX, KSWY, K, L, 
     *    MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4,
     *    AIT, BIT, CIT, DIT, MIT, NIT, IS


      integer I, J

      Logical SINGLR

      real*8 ALPHA
      real*8 BETA
      real*8 GAMA
      real*8 XNU



      real*8 YJ, DJ, EJ, FJ
      real*8 XI, AI, BI, CI

      SINGLR = .FALSE.
c
c     check if the boundary conditions are
c     entirely periodic and/or mixed
c
      If ((MBDCND.NE.0.AND.MBDCND.NE.3).OR.(NBDCND.NE.0.AND.NBDCND.NE.3)
     *    ) Return
c
c     check that mixed conditions are pure neuman
c
      If (MBDCND.EQ.3) Then
        If (ALPHA.NE.0.0.OR.BETA.NE.0.0) Return
      End If
      If (NBDCND.EQ.3) Then
        If (GAMA.NE.0.0.OR.XNU.NE.0.0) Return
      End If
c
c     check that non-derivative coefficient functions
c     are zero
c
      Do 10 I = IS, MS
        XI = AIT + DBLE(I-1) * DLX
        Call COFX(XI,AI,BI,CI)
        If (CI.NE.0.0) Return
   10 Continue
      Do 20 J = JS, NS
        YJ = CIT + DBLE(J-1) * DLY
        Call COFY(YJ,DJ,EJ,FJ)
        If (FJ.NE.0.0) Return
   20 Continue
c
c     the operator must be singular if this point is reached
c
      SINGLR = .TRUE.
      Return
      End


















c      Subroutine DEFER(COFX,COFY,IDMN,USOL,GRHS)
      Subroutine DEFER(IDMN,USOL,GRHS)
c
c     this subroutine first approximates the truncation error given by
c     trun1(x,y)=dlx**2*tx+dly**2*ty where
c     tx=afun(x)*uxxxx/12.0+bfun(x)*uxxx/6.0 on the interior and
c     at the boundaries if periodic(here uxxx,uxxxx are the third
c     and fourth partial derivatives of u with respect to x).
c     tx is of the form afun(x)/3.0*(uxxxx/4.0+uxxx/dlx)
c     at x=a or x=b if the boundary condition there is mixed.
c     tx=0.0 along specified boundaries.  ty has symmetric form
c     in y with x,afun(x),bfun(x) replaced by y,dfun(y),efun(y).
c     the second order solution in usol is used to approximate
c     (via second order finite differencing) the truncation error
c     and the result is added to the right hand side in grhs
c     and then transferred to usol to be used as a new right
c     hand side when calling blktri for a fourth order solution.
c
      IMPLICIT NONE
      integer KSWX, KSWY, K 
      integer AIT, BIT, CIT, DIT, IS, MIT, NIT
      real*8 MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4

      integer L
      Common /SPLP/ KSWX, KSWY, K, L, 
     *    MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4,
     *    AIT, BIT, CIT, DIT, MIT, NIT, IS


      integer IDMN
      Dimension GRHS(IDMN,1), USOL(IDMN,1)
      real*8 GRHS, USOL, UXXX, UXXXX, UYYY, UYYYY
c      External COFX, COFY

      real*8 YJ, DJ, EJ, FJ
      real*8 XI, AI, BI, CI, TX, TY

      integer I, J
c
c     compute truncation error approximation over the entire mesh
c
      Do 20 J = JS, NS
        YJ = CIT + DBLE(J-1) * DLY
        Call COFY(YJ,DJ,EJ,FJ)
        Do 10 I = IS, MS
          XI = AIT + DBLE(I-1) * DLX
          Call COFX(XI,AI,BI,CI)
c
c     compute partial derivative approximations at (xi,yj)
c
          Call SEPDX(USOL,IDMN,I,J,UXXX,UXXXX)
          Call SEPDY(USOL,IDMN,I,J,UYYY,UYYYY)
          TX = AI * UXXXX / 12.0 + BI * UXXX / 6.0
          TY = DJ * UYYYY / 12.0 + EJ * UYYY / 6.0
c
c     reset form of truncation if at boundary which is non-periodic
c
          If (.NOT.(KSWX.EQ.1.OR.(I.GT.1.AND.I.LT.K))) Then
            TX = AI / 3.0 * (UXXXX/4.0+UXXX/DLX)
          End If
          If (.NOT.(KSWY.EQ.1.OR.(J.GT.1.AND.J.LT.L))) Then
            TY = DJ / 3.0 * (UYYYY/4.0+UYYY/DLY)
          End If
          GRHS(I,J) = GRHS(I,J) + DLX ** 2 * TX + DLY ** 2 * TY
   10   Continue
   20 Continue
c
c     reset the right hand side in usol
c
      Do 40 I = IS, MS
        Do 30 J = JS, NS
          USOL(I,J) = GRHS(I,J)
   30   Continue
   40 Continue
      Return
c
c revision history---
c
c december 1979    first added to nssl
c-----------------------------------------------------------------------
      End











      Subroutine BLKTRI(IFLG,NP,N,AN,BN,CN,MP,M,AM,BM,CM,IDIMY,Y,IERROR,
     *    W)
c
c dimension of           an(n),bn(n),cn(n),am(m),bm(m),cm(m),y(idimy,n),
c arguments              w(see argument list)
c
c latest revision        january 1985
c
c usage                  call blktri (iflg,np,n,an,bn,cn,mp,m,am,bm,
c                                     cm,idimy,y,ierror,w)
c
c purpose                blktri solves a system of linear equations
c                        of the form
c
c                        an(j)*x(i,j-1) + am(i)*x(i-1,j) +
c                        (bn(j)+bm(i))*x(i,j) + cn(j)*x(i,j+1) +
c                        cm(i)*x(i+1,j) = y(i,j)
c
c                        for i = 1,2,...,m  and  j = 1,2,...,n.
c
c                        i+1 and i-1 are evaluated modulo m and
c                        j+1 and j-1 modulo n, i.e.,
c
c                        x(i,0) = x(i,n),  x(i,n+1) = x(i,1),
c                        x(0,j) = x(m,j),  x(m+1,j) = x(1,j).
c
c                        these equations usually result from the
c                        discretization of separable elliptic
c                        equations.  boundary conditions may be
c                        dirichlet, neumann, or periodic.
c
c arguments
c
c on input               iflg
c
c                          = 0  initialization only.
c                               certain quantities that depend on np,
c                               n, an, bn, and cn are computed and
c                               stored in the work array w.
c
c                          = 1  the quantities that were computed
c                               in the initialization are used
c                               to obtain the solution x(i,j).
c
c                               note:
c                               a call with iflg=0 takes
c                               approximately one half the time
c                               as a call with iflg = 1.
c                               however, the initialization does
c                               not have to be repeated unless np,
c                               n, an, bn, or cn change.
c
c                        np
c                          = 0  if an(1) and cn(n) are not zero,
c                               which corresponds to periodic
c                               bounary conditions.
c
c                          = 1  if an(1) and cn(n) are zero.
c
c                        n
c                          the number of unknowns in the j-direction.
c                          n must be greater than 4.
c                          the operation count is proportional to
c                          mnlog2(n), hence n should be selected
c                          less than or equal to m.
c
c                        an,bn,cn
c                          one-dimensional arrays of length n
c                          that specify the coefficients in the
c                          linear equations given above.
c
c                        mp
c                          = 0  if am(1) and cm(m) are not zero,
c                               which corresponds to periodic
c                               boundary conditions.
c
c                          = 1  if am(1) = cm(m) = 0  .
c
c                        m
c                          the number of unknowns in the i-direction.
c                           m must be greater than 4.
c
c                        am,bm,cm
c                          one-dimensional arrays of length m that
c                          specify the coefficients in the linear
c                          equations given above.
c
c                        idimy
c                          the row (or first) dimension of the
c                          two-dimensional array y as it appears
c                          in the program calling blktri.
c                          this parameter is used to specify the
c                          variable dimension of y.
c                          idimy must be at least m.
c
c                        y
c                          a two-dimensional array that specifies
c                          the values of the right side of the linear
c                          system of equations given above.
c                          y must be dimensioned at least m*n.
c
c                        w
c                          a one-dimensional array that must be
c                          provided by the user for work space.
c                          if np=1 define k=int(log2(n))+1 and
c                          set l=2**(k+1) then w must have dimension
c                          (k-2)*l+k+5+max(2n,6m)
c
c                          if np=0 define k=int(log2(n-1))+1 and
c                          set l=2**(k+1) then w must have dimension
c                          (k-2)*l+k+5+2n+max(2n,6m)
c
c                          **important**
c                          for purposes of checking, the required
c                          dimension of w is computed by blktri and
c                          stored in w(1) in floating point format.
c
c arguments
c
c on output              y
c                          contains the solution x.
c
c                        ierror
c                          an error flag that indicates invalid
c                          input parameters.  except for number zer0,
c                          a solution is not attempted.
c
c                        = 0  no error.
c                        = 1  m is less than 5
c                        = 2  n is less than 5
c                        = 3  idimy is less than m.
c                        = 4  blktri failed while computing results
c                             that depend on the coefficient arrays
c                             an, bn, cn.  check these arrays.
c                        = 5  an(j)*cn(j-1) is less than 0 for some j.
c
c                             possible reasons for this condition are
c                             1. the arrays an and cn are not correct
c                             2. too large a grid spacing was used
c                                in the discretization of the elliptic
c                                equation.
c                             3. the linear equations resulted from a
c                                partial differential equation which
c                                was not elliptic.
c
c                        w
c                           contains intermediate values that must
c                           not be destroyed if blktri will be called
c                           again with iflg=1. w(1) contains the
c                           number of locations required by w in
c                           floating point format.
c
c
c special conditions     the algorithm may fail if abs(bm(i)+bn(j))
c                        is less than abs(am(i))+abs(an(j))+
c                        abs(cm(i))+abs(cn(j))
c                        for some i and j. the algorithm will also
c                        fail if an(j)*cn(j-1) is less than zero for
c                        some j.
c                        see the description of the output parameter
c                        ierror.
c
c i/o                    none
c
c precision              single
c
c required library       comf and q8qst4, which are automatically loaded
c files                  ncar's cray machines.
c
c language               fortran
c
c history                written by paul swarztrauber at ncar in the
c                        early 1970's.  rewritten and released in
c                        january, 1980.
c
c algorithm              generalized cyclic reduction
c
c portability            fortran 66.  approximate machine accuracy
c                        is computed in function epmach.
c
c references             swarztrauber,p. and r. sweet, 'efficient
c                        fortran subprograms for the solution of
c                        elliptic equations'
c                        ncar tn/ia-109, july, 1975, 138 pp.
c
c                        swarztrauber p. n.,a direct method for
c                        the discrete solution of separable
c                        elliptic equations, s.i.a.m.
c                        j. numer. anal.,11(1974) pp. 1136-1150.
c***********************************************************************

      IMPLICIT NONE

      integer IDIMY
      Dimension AN(1), BN(1), CN(1), AM(1), BM(1), CM(1), Y(IDIMY,1) 
      Dimension W(2)
      External PROD, PRODP, CPROD, CPRODP

      integer M, N, NP, MP
      integer K, NM, NCMPLX, IK, NL
      real*8 EPS, NPP, CNV
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS

      real*8 AN, BN, CN, AM, BM, CM, Y, W



      integer IERROR, NH, IWAH, IW1, IWBH, IW2, IW3, IWD, IWW, IWU
      integer IFLG

c
c the following call is for monitoring library use at ncar
c
      Logical Q8Q4
      Save Q8Q4
      Data Q8Q4 /.TRUE./
      If (Q8Q4) Then
c         call q8qst4('loclib','blktri','blktri','version 01')
        Q8Q4 = .FALSE.
      End If
c
c test m and n for the proper form
c
      NM = N
      IERROR = 0
      If (M.LT.5) Then
        IERROR = 1
      Else
        If (NM.LT.3) Then
          IERROR = 2
        Else
          If (IDIMY.LT.M) Then
            IERROR = 3
          Else
            NH = N
            NPP = NP
            If (NPP.NE.0) Then
              NH = NH + 1
            End If
            IK = 2
            K = 1
   10       IK = IK + IK
            K = K + 1
            If (NH.GT.IK) Go To 10
            NL = IK
            IK = IK + IK
            NL = NL - 1
            IWAH = (K-2) * IK + K + 6
            If (NPP.NE.0) Then
c
c     divide w into working sub arrays
c
              IW1 = IWAH
              IWBH = IW1 + NM
              W(1) = DBLE(IW1-1+DMAX1(2.d0*NM,6.d0*M))
            Else
              IWBH = IWAH + NM + NM
              IW1 = IWBH
              W(1) = DBLE(IW1-1+DMAX1(2.d0*NM,6.d0*M))
              NM = NM - 1
            End If
c
c subroutine comp b computes the roots of the b polynomials
c
            If (IERROR.EQ.0) Then
              IW2 = IW1 + M
              IW3 = IW2 + M
              IWD = IW3 + M
              IWW = IWD + M
              IWU = IWW + M
              If (IFLG.EQ.0) Then
                Call COMPB(NL,IERROR,AN,BN,CN,W(2),W(IWAH),W(IWBH))
              Else
                If (MP.NE.0) Then
c
c subroutine blktr1 solves the linear system
c
                  Call BLKTR1(NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),
     *                W(IW1),W(IW2),W(IW3),W(IWD),W(IWW),W(IWU),1)
c                  Call BLKTR1(NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),
c     *                W(IW1),W(IW2),W(IW3),W(IWD),W(IWW),W(IWU),PROD,
c     *                CPROD)
                Else
c                  Call BLKTR1(NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),
c     *                W(IW1),W(IW2),W(IW3),W(IWD),W(IWW),W(IWU),PRODP,
c     *                CPRODP)
                  Call BLKTR1(NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),
     *                W(IW1),W(IW2),W(IW3),W(IWD),W(IWW),W(IWU),2)
                End If
              End If
            End If
          End If
        End If
      End If
      Return
      End









      Subroutine BLKTR1(N,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,B,W1,W2,W3,WD,WW,
c         WU,PRDCT,CPRDCT)
     *    WU, PFLAG )
c
c blktr1 solves the linear system
c
c b  contains the roots of all the b polynomials
c w1,w2,w3,wd,ww,wu  are all working arrays
c prdct  is either prodp or prod depending on whether the boundary
c conditions in the m direction are periodic or not
c cprdct is either cprodp or cprod which are the complex versions
c of prodp and prod. these are called in the event that some
c of the roots of the b sub p polynomial are complex
c
c If PFLAG==1, then use PROD and CPROD as callable functions. 
c Otherwise use PRODP and CPRODP
c

      IMPLICIT NONE


      integer IDIMY
      Dimension AN(1), BN(1), CN(1), AM(1), BM(1), CM(1), B(1), W1(1), 
     *    W2(1), W3(1), WD(1), WW(1), WU(1), Y(IDIMY,1)
      real*8 AN, BN, CN, AM, BM, CM, B, W1
      real*8 W2, W3, WD, WW, WU, Y
      real*8 DUM
      integer M, N
      integer PFLAG


      integer K, NM, NCMPLX, IK
      real*8 EPS, NPP, CNV
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS


      integer  I2,IR,IM2,NM2, I1, IRM1, IM3, NM3
      integer I3, IM1, NM1, KDO, L, I4, IF,I, IPI1, IPI2, IPI3
      integer NC, IDXA, IP2, NA, NP2, IP1, NP1, J, IP3, NP3
      integer IZ, NZ, IDXC, IZR, LL, IFD, IP, NP, IMI1, IMI2

c
c Added this to avoid compiler warning on g77
c John Evans
      Complex cbip, cw1, cw3, cww


c
c begin reduction phase
c
      KDO = K - 1
      Do 40 L = 1, KDO
        IR = L - 1
        I2 = 2 ** IR
        I1 = I2 / 2
        I3 = I2 + I1
        I4 = I2 + I2
        IRM1 = IR - 1
        Call INDXB(I2,IR,IM2,NM2)
        Call INDXB(I1,IRM1,IM3,NM3)
        Call INDXB(I3,IRM1,IM1,NM1)

        If (PFLAG.eq.1) Then
            Call PROD(NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,Y(1,I2),
     *                 W3, M, AM,BM,CM,WD,WW,WU)
        Else
            Call PRODP(NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,Y(1,I2),
     *                 W3,M,AM,BM,CM,WD,WW,WU)
        End If
c
c        Call PRDCT(NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,Y(1,I2),W3,M,
c     *      AM,BM,CM,WD,WW,WU)
c

        IF = 2 ** K
        Do 30 I = I4, IF, I4
          If (I.LE.NM) Then
            IPI1 = I + I1
            IPI2 = I + I2
            IPI3 = I + I3
            Call INDXC(I,IR,IDXC,NC)
            If (I.LT.IF) Then
              Call INDXA(I,IR,IDXA,NA)
              Call INDXB(I-I1,IRM1,IM1,NM1)
              Call INDXB(IPI2,IR,IP2,NP2)
              Call INDXB(IPI1,IRM1,IP1,NP1)
              Call INDXB(IPI3,IRM1,IP3,NP3)

c              Call PRDCT(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W3,W1,M,AM,
c     *            BM,CM,WD,WW,WU)

              If (PFLAG .eq.1) Then
                Call PROD(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W3,W1,
     *                   M,AM,BM,CM,WD,WW,WU)
              Else
                Call PRODP(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W3,W1,
     *                    M,AM,BM,CM,WD,WW,WU)
              End If


              If (IPI2.GT.NM) Then
                Do 10 J = 1, M
                  W3(J) = 0.
                  W2(J) = 0.
   10           Continue
              Else

c                Call PRDCT(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,
c     *              Y(1,IPI2),W3,M,AM,BM,CM,WD,WW,WU)
c                Call PRDCT(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W3,W2,M,
c     *              AM,BM,CM,WD,WW,WU)

                If (PFLAG.eq.1) Then
                    Call PROD(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,
     *                  Y(1,IPI2),W3,M,AM,BM,CM,WD,WW,WU)
                    Call PROD(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W3,
     *                   W2,M,AM,BM,CM,WD,WW,WU)
                Else
                    Call PRODP(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,
     *                  Y(1,IPI2),W3,M,AM,BM,CM,WD,WW,WU)
                    Call PRODP(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W3,
     *                   W2,M,AM,BM,CM,WD,WW,WU)
                End If

              End If
              Do 20 J = 1, M
                Y(J,I) = W1(J) + W2(J) + Y(J,I)
   20         Continue
            End If
          End If
   30   Continue
   40 Continue
      If (NPP.EQ.0) Then
c
c     the periodic case is treated using the capacitance matrix method
c
        IF = 2 ** K
        I = IF / 2
        I1 = I / 2
        Call INDXB(I-I1,K-2,IM1,NM1)
        Call INDXB(I+I1,K-2,IP1,NP1)
        Call INDXB(I,K-1,IZ,NZ)



c        Call PRDCT(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,Y(1,I),W1,M,AM,
c     *      BM,CM,WD,WW,WU)
        If (PFLAG.eq.1) Then
            Call PROD(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,Y(1,I),W1,
     *                M,AM,BM,CM,WD,WW,WU)
        Else
            Call PRODP(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,Y(1,I),
     *                W1,M,AM,BM,CM,WD,WW,WU)
        End If


        IZR = I
        Do 50 J = 1, M
          W2(J) = W1(J)
   50   Continue
        Do 70 LL = 2, K
          L = K - LL + 1
          IR = L - 1
          I2 = 2 ** IR
          I1 = I2 / 2
          I = I2
          Call INDXC(I,IR,IDXC,NC)
          Call INDXB(I,IR,IZ,NZ)
          Call INDXB(I-I1,IR-1,IM1,NM1)
          Call INDXB(I+I1,IR-1,IP1,NP1)



c          Call PRDCT(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W1,W1,M,AM,BM,
c     *        CM,WD,WW,WU)
          If (PFLAG.eq.1) Then
              Call PROD(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W1,W1,M,
     *                  AM,BM,CM,WD,WW,WU)
          Else
              Call PRODP(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W1,W1,
     *                  M,AM,BM,CM,WD,WW,WU)
          End If



          Do 60 J = 1, M
            W1(J) = Y(J,I) + W1(J)
   60     Continue



c          Call PRDCT(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,W1,M,AM,BM,
c     *        CM,WD,WW,WU)
          If (PFLAG.eq.1) Then
              Call PROD(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,W1,M,
     *                  AM,BM,CM,WD,WW,WU)
          Else
              Call PRODP(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,W1,
     *                  M,AM,BM,CM,WD,WW,WU)
          End If



   70   Continue
        Do 100 LL = 2, K
          L = K - LL + 1
          IR = L - 1
          I2 = 2 ** IR
          I1 = I2 / 2
          I4 = I2 + I2
          IFD = IF - I2
          Do 90 I = I2, IFD, I4
            If (I-I2-IZR.EQ.0) Then
              If (I.GT.NM) Go To 100
              Call INDXA(I,IR,IDXA,NA)
              Call INDXB(I,IR,IZ,NZ)
              Call INDXB(I-I1,IR-1,IM1,NM1)
              Call INDXB(I+I1,IR-1,IP1,NP1)



c              Call PRDCT(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W2,W2,M,AM,
c     *            BM,CM,WD,WW,WU)
              If (PFLAG.eq.1) Then
                  Call PROD(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W2,W2,
     *                      M,AM,BM,CM,WD,WW,WU)
              Else
                  Call PRODP(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W2,
     *                      W2,M,AM,BM,CM,WD,WW,WU)
              End If



              Do 80 J = 1, M
                W2(J) = Y(J,I) + W2(J)
   80         Continue




c              Call PRDCT(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W2,W2,M,
c     *            AM,BM,CM,WD,WW,WU)
              If (PFLAG.eq.1) Then
                  Call PROD(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W2,
     *                      W2,M,AM,BM,CM,WD,WW,WU)
              Else
                  Call PRODP(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,
     *                       W2,W2,M,AM,BM,CM,WD,WW,WU)
              End If




              IZR = I
              If (I.EQ.NM) Go To 110
            End If
   90     Continue
  100   Continue
  110   Do 120 J = 1, M
          Y(J,NM+1) = Y(J,NM+1) - CN(NM+1) * W1(J) - AN(NM+1) * W2(J)
  120   Continue
        Call INDXB(IF/2,K-1,IM1,NM1)
        Call INDXB(IF,K-1,IP,NP)



        If (NCMPLX.NE.0) Then

          cbip = CMPLX ( B(IP) )
          cw1 = CMPLX ( W1(1) )
          cw3 = CMPLX ( W3(1) )
          cww = CMPLX ( WW(1) )

c          Call CPRDCT(NM+1,B(IP),NM1,B(IM1),0,DUM,0,DUM,Y(1,NM+1),
c     *        Y(1,NM+1),M,AM,BM,CM,W1,W3,WW)
          If (PFLAG.eq.1) Then
              Call CPROD(NM+1, cbip, NM1,B(IM1),0,DUM,0,DUM,
     *                   Y(1,NM+1),Y(1,NM+1),M,AM,BM,CM, cw1, cw3, cww)
          Else
              Call CPRODP(NM+1,cbip,NM1,B(IM1),0,DUM,0,DUM,
     *                    Y(1,NM+1),Y(1,NM+1),M,AM,BM,CM, cw1, cw3, cww)
          End If


        Else

c          Call PRDCT(NM+1,B(IP),NM1,B(IM1),0,DUM,0,DUM,Y(1,NM+1),
c     *        Y(1,NM+1),M,AM,BM,CM,WD,WW,WU)
          If (PFLAG.eq.1) Then
              Call PROD(NM+1,B(IP),NM1,B(IM1),0,DUM,0,DUM,Y(1,NM+1),
     *            Y(1,NM+1),M,AM,BM,CM,WD,WW,WU)
          Else
              Call PRODP(NM+1,B(IP),NM1,B(IM1),0,DUM,0,DUM,Y(1,NM+1),
     *            Y(1,NM+1),M,AM,BM,CM,WD,WW,WU)
          End If


        End If




        Do 130 J = 1, M
          W1(J) = AN(1) * Y(J,NM+1)
          W2(J) = CN(NM) * Y(J,NM+1)
          Y(J,1) = Y(J,1) - W1(J)
          Y(J,NM) = Y(J,NM) - W2(J)
  130   Continue
        Do 150 L = 1, KDO
          IR = L - 1
          I2 = 2 ** IR
          I4 = I2 + I2
          I1 = I2 / 2
          I = I4
          Call INDXA(I,IR,IDXA,NA)
          Call INDXB(I-I2,IR,IM2,NM2)
          Call INDXB(I-I2-I1,IR-1,IM3,NM3)
          Call INDXB(I-I1,IR-1,IM1,NM1)


c         Call PRDCT(NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,W1,W1,M,AM,
c    *        BM,CM,WD,WW,WU)
c         Call PRDCT(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W1,W1,M,AM,BM,
c    *        CM,WD,WW,WU)
          If (PFLAG.eq.1) Then
              Call PROD(NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,W1,W1,
     *                  M,AM,BM,CM,WD,WW,WU)
              Call PROD(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W1,W1,M,
     *                  AM,BM,CM,WD,WW,WU)
          Else
              Call PRODP(NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,W1,
     *                   W1,M,AM,BM,CM,WD,WW,WU)
              Call PRODP(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W1,W1,
     *                   M,AM,BM,CM,WD,WW,WU)
          End If




          Do 140 J = 1, M
            Y(J,I) = Y(J,I) - W1(J)
  140     Continue
  150   Continue
c
        IZR = NM
        Do 180 L = 1, KDO
          IR = L - 1
          I2 = 2 ** IR
          I1 = I2 / 2
          I3 = I2 + I1
          I4 = I2 + I2
          IRM1 = IR - 1
          Do 170 I = I4, IF, I4
            IPI1 = I + I1
            IPI2 = I + I2
            IPI3 = I + I3
            If (IPI2.NE.IZR) Then
              If (I-IZR) 170, 180, 170
            End If
            Call INDXC(I,IR,IDXC,NC)
            Call INDXB(IPI2,IR,IP2,NP2)
            Call INDXB(IPI1,IRM1,IP1,NP1)
            Call INDXB(IPI3,IRM1,IP3,NP3)




c            Call PRDCT(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,W2,W2,M,
c     *          AM,BM,CM,WD,WW,WU)
c            Call PRDCT(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W2,W2,M,AM,BM,
c     *          CM,WD,WW,WU)
            If (PFLAG.eq.1) Then
                Call PROD(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,W2,
     *                    W2,M,AM,BM,CM,WD,WW,WU)
                Call PROD(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W2,W2,M,
     *                    AM,BM,CM,WD,WW,WU)
            Else
                Call PRODP(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,W2,
     *                     W2,M,AM,BM,CM,WD,WW,WU)
                Call PRODP(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W2,W2,M,
     *                     AM,BM,CM,WD,WW,WU)
            End If



            Do 160 J = 1, M
              Y(J,I) = Y(J,I) - W2(J)
  160       Continue
            IZR = I
            Go To 180
  170     Continue
  180   Continue
      End If
c
c begin back substitution phase
c
      Do 230 LL = 1, K
        L = K - LL + 1
        IR = L - 1
        IRM1 = IR - 1
        I2 = 2 ** IR
        I1 = I2 / 2
        I4 = I2 + I2
        IFD = IF - I2
        Do 220 I = I2, IFD, I4
          If (I.LE.NM) Then
            IMI1 = I - I1
            IMI2 = I - I2
            IPI1 = I + I1
            IPI2 = I + I2
            Call INDXA(I,IR,IDXA,NA)
            Call INDXC(I,IR,IDXC,NC)
            Call INDXB(I,IR,IZ,NZ)
            Call INDXB(IMI1,IRM1,IM1,NM1)
            Call INDXB(IPI1,IRM1,IP1,NP1)



            If (I.LE.I2) Then
              Do 190 J = 1, M
                W1(J) = 0.
  190         Continue

            Else

c              Call PRDCT(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),Y(1,IMI2),
c     *            W1,M,AM,BM,CM,WD,WW,WU)
              If (PFLAG.eq.1) Then
                  Call PROD(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),
     *                      Y(1,IMI2),W1,M,AM,BM,CM,WD,WW,WU)
              Else
                  Call PRODP(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),
     *                       Y(1,IMI2),W1,M,AM,BM,CM,WD,WW,WU)
              End If

            End If



            If (IPI2.GT.NM) Then

              Do 200 J = 1, M
                W2(J) = 0.
  200         Continue

            Else

c              Call PRDCT(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),Y(1,IPI2),
c     *            W2,M,AM,BM,CM,WD,WW,WU)
              If (PFLAG.eq.1) Then
                  Call PROD(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),
     *                      Y(1,IPI2),W2,M,AM,BM,CM,WD,WW,WU)
              Else
                  Call PRODP(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),
     *                       Y(1,IPI2),W2,M,AM,BM,CM,WD,WW,WU)
              End If


            End If



            Do 210 J = 1, M
              W1(J) = Y(J,I) + W1(J) + W2(J)
  210       Continue


c            Call PRDCT(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,Y(1,I),M,
c     *          AM,BM,CM,WD,WW,WU)
            If (PFLAG.eq.1) Then
                Call PROD(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,
     *                    Y(1,I),M,AM,BM,CM,WD,WW,WU)
            Else
                Call PRODP(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,
     *                     Y(1,I),M,AM,BM,CM,WD,WW,WU)
            End If


          End If


  220   Continue
  230 Continue
      Return
      End

















      Subroutine INDXB(I,IR,IDX,IDP)
c
c b(idx) is the location of the first root of the b(i,ir) polynomial
c
      IMPLICIT NONE

      integer IDP, IDX, I, IR
      integer IZH, ID, IPL

      integer K, NM, NCMPLX, IK
      real*8 EPS, NPP, CNV
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS
      IDP = 0
      If (IR) 30, 10, 20
   10 If (I.GT.NM) Go To 30
      IDX = I
      IDP = 1
      Return
   20 IZH = 2 ** IR
      ID = I - IZH - IZH
      IDX = ID + ID + (IR-1) * IK + IR + (IK-I) / IZH + 4
      IPL = IZH - 1
      IDP = IZH + IZH - 1
      If (I-IPL-NM.GT.0) Then
        IDP = 0
        Return
      End If
      If (I+IPL-NM.GT.0) Then
        IDP = NM + IPL - I + 1
      End If
   30 Return
      End








      Subroutine INDXA(I,IR,IDXA,NA)
      IMPLICIT NONE

      integer K, NM, NCMPLX, IK
      real*8 EPS, NPP, CNV
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS

      integer NA, IDXA, I, IR



      NA = 2 ** IR
      IDXA = I - NA + 1
      If (I.GT.NM) Then
        NA = 0
      End If
      Return
      End





      Subroutine INDXC(I,IR,IDXC,NC)

      IMPLICIT NONE
      integer K, NM, NCMPLX, IK, NC, IR, I, IDXC
      real*8 EPS, NPP, CNV
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS

      NC = 2 ** IR
      IDXC = I
      If (IDXC+NC-1-NM.GT.0) Then
        NC = 0
      End If
      Return
      End





      Subroutine PROD(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,W,U)
c
c prod applies a sequence of matrix operations to the vector x and
c stores the result in y
c bd,bm1,bm2 are arrays containing roots of certian b polynomials
c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
c aa   array containing scalar multipliers of the vector x
c na is the length of the array aa
c x,y  the matrix operations are applied to x and the result is y
c a,b,c  are arrays which contain the tridiagonal matrix
c m  is the order of the matrix
c d,w,u are working arrays
c is  determines whether or not a change in sign is made
c
      IMPLICIT NONE

      Dimension A(1), B(1), C(1), X(1), Y(1), D(2), W(2), BD(1), 
     *    BM1(1), BM2(1), AA(1), U(1)

      real*8 bd,bm1, bm2, aa, a, b, c, d, u, w, x, y, rt
      real*8 den
      integer nd, nm1, nm2, na, m, j, mm, id, ibr, m1, m2, ia
      integer k



      Do 10 J = 1, M
        W(J) = X(J)
        Y(J) = W(J)
   10 Continue
      MM = M - 1
      ID = ND
      IBR = 0
      M1 = NM1
      M2 = NM2
      IA = NA
   20 If (IA.GT.0) Then
        RT = AA(IA)
        If (ND.EQ.0) RT = -RT
        IA = IA - 1
c
c scalar multiplication
c
        Do 30 J = 1, M
          Y(J) = RT * W(J)
   30   Continue
      End If
      If (ID.GT.0) Then
        RT = BD(ID)
        ID = ID - 1
        If (ID.EQ.0) IBR = 1
c
c begin solution to system
c
        D(M) = A(M) / (B(M)-RT)
        W(M) = Y(M) / (B(M)-RT)
        Do 40 J = 2, MM
          K = M - J
          DEN = B(K+1) - RT - C(K+1) * D(K+2)
          D(K+1) = A(K+1) / DEN
          W(K+1) = (Y(K+1)-C(K+1)*W(K+2)) / DEN
   40   Continue
        DEN = B(1) - RT - C(1) * D(2)
        W(1) = 1.
        If (DEN.NE.0) Then
          W(1) = (Y(1)-C(1)*W(2)) / DEN
        End If
        Do 50 J = 2, M
          W(J) = W(J) - D(J) * W(J-1)
   50   Continue
        If (NA) 80, 80, 20
   60   Do 70 J = 1, M
          Y(J) = W(J)
   70   Continue
        IBR = 1
        Go To 20
   80   If (M1.LE.0) Then
          If (M2) 60, 60, 90
        End If
        If (M2.GT.0) Then
          If (ABS(BM1(M1)).LE.ABS(BM2(M2))) Go To 90
        End If
        If (IBR.LE.0) Then
          If (ABS(BM1(M1)-BD(ID)).LT.ABS(BM1(M1)-RT)) Go To 60
        End If
        RT = RT - BM1(M1)
        M1 = M1 - 1
        Go To 100
   90   If (IBR.LE.0) Then
          If (ABS(BM2(M2)-BD(ID)).LT.ABS(BM2(M2)-RT)) Go To 60
        End If
        RT = RT - BM2(M2)
        M2 = M2 - 1
  100   Do 110 J = 1, M
          Y(J) = Y(J) + RT * W(J)
  110   Continue
        Go To 20
      End If
      Return
      End





      Subroutine PRODP(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,U,W)
c
c prodp applies a sequence of matrix operations to the vector x and
c stores the result in y        periodic boundary conditions
c
c bd,bm1,bm2 are arrays containing roots of certian b polynomials
c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
c aa   array containing scalar multipliers of the vector x
c na is the length of the array aa
c x,y  the matrix operations are applied to x and the result is y
c a,b,c  are arrays which contain the tridiagonal matrix
c m  is the order of the matrix
c d,u,w are working arrays
c is  determines whether or not a change in sign is made
c
      IMPLICIT NONE



      Dimension A(1), B(1), C(1), X(1), Y(1), D(1), U(1), BD(1), 
     *    BM1(1), BM2(1), AA(1), W(1)


      real*8 bd,bm1, bm2, aa, a, b, c, d, u, w, x, y, rt, bh, ym, v
      real*8 den, am
      integer nd, nm1, nm2, na, m, j, mm, mm2, id, ibr, m1, m2, ia
      integer k


      Do 10 J = 1, M
        Y(J) = X(J)
        W(J) = Y(J)
   10 Continue
      MM = M - 1
      MM2 = M - 2
      ID = ND
      IBR = 0
      M1 = NM1
      M2 = NM2
      IA = NA
   20 If (IA.GT.0) Then
        RT = AA(IA)
        If (ND.EQ.0) RT = -RT
        IA = IA - 1
        Do 30 J = 1, M
          Y(J) = RT * W(J)
   30   Continue
      End If
      If (ID.GT.0) Then
        RT = BD(ID)
        ID = ID - 1
        If (ID.EQ.0) IBR = 1
c
c begin solution to system
c
        BH = B(M) - RT
        YM = Y(M)
        DEN = B(1) - RT
        D(1) = C(1) / DEN
        U(1) = A(1) / DEN
        W(1) = Y(1) / DEN
        V = C(M)
        If (MM2.GE.2) Then
          Do 40 J = 2, MM2
            DEN = B(J) - RT - A(J) * D(J-1)
            D(J) = C(J) / DEN
            U(J) = -A(J) * U(J-1) / DEN
            W(J) = (Y(J)-A(J)*W(J-1)) / DEN
            BH = BH - V * U(J-1)
            YM = YM - V * W(J-1)
            V = -V * D(J-1)
   40     Continue
        End If
        DEN = B(M-1) - RT - A(M-1) * D(M-2)
        D(M-1) = (C(M-1)-A(M-1)*U(M-2)) / DEN
        W(M-1) = (Y(M-1)-A(M-1)*W(M-2)) / DEN
        AM = A(M) - V * D(M-2)
        BH = BH - V * U(M-2)
        YM = YM - V * W(M-2)
        DEN = BH - AM * D(M-1)
        If (DEN.NE.0) Then
          W(M) = (YM-AM*W(M-1)) / DEN
        Else
          W(M) = 1.
        End If
        W(M-1) = W(M-1) - D(M-1) * W(M)
        Do 50 J = 2, MM
          K = M - J
          W(K) = W(K) - D(K) * W(K+1) - U(K) * W(M)
   50   Continue
        If (NA) 80, 80, 20
   60   Do 70 J = 1, M
          Y(J) = W(J)
   70   Continue
        IBR = 1
        Go To 20
   80   If (M1.LE.0) Then
          If (M2) 60, 60, 90
        End If
        If (M2.GT.0) Then
          If (ABS(BM1(M1)).LE.ABS(BM2(M2))) Go To 90
        End If
        If (IBR.LE.0) Then
          If (ABS(BM1(M1)-BD(ID)).LT.ABS(BM1(M1)-RT)) Go To 60
        End If
        RT = RT - BM1(M1)
        M1 = M1 - 1
        Go To 100
   90   If (IBR.LE.0) Then
          If (ABS(BM2(M2)-BD(ID)).LT.ABS(BM2(M2)-RT)) Go To 60
        End If
        RT = RT - BM2(M2)
        M2 = M2 - 1
  100   Do 110 J = 1, M
          Y(J) = Y(J) + RT * W(J)
  110   Continue
        Go To 20
      End If
      Return
      End

















      Subroutine CPROD(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,YY,M,A,B,C,D,W,Y)
c
c prod applies a sequence of matrix operations to the vector x and
c stores the result in yy           (complex case)
c aa   array containing scalar multipliers of the vector x
c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
c bd,bm1,bm2 are arrays containing roots of certian b polynomials
c na is the length of the array aa
c x,yy the matrix operations are applied to x and the result is yy
c a,b,c  are arrays which contain the tridiagonal matrix
c m  is the order of the matrix
c d,w,y are working arrays
c isgn  determines whether or not a change in sign is made
c
      IMPLICIT NONE
      Complex Y, D, W, BD, CRT, DEN, Y1, Y2
      Dimension A(1), B(1), C(1), X(1), Y(2), D(2), W(2), BD(1), 
     *    BM1(1), BM2(1), AA(1), YY(1)



      integer J, M, MM, ND, NM1, NM2, ID, M1, M2, IA
      integer NA, IFLG, K
      real*8 BM1, BM2, AA, X, YY, A, B, C, RT


      Do 10 J = 1, M
        Y(J) = DCMPLX(X(J),0.d0)
   10 Continue
      MM = M - 1
      ID = ND
      M1 = NM1
      M2 = NM2
      IA = NA
   20 IFLG = 0
      If (ID.GT.0) Then
        CRT = BD(ID)
        ID = ID - 1
c
c begin solution to system
c
        D(M) = A(M) / (B(M)-CRT)
        W(M) = Y(M) / (B(M)-CRT)
        Do 30 J = 2, MM
          K = M - J
          DEN = B(K+1) - CRT - C(K+1) * D(K+2)
          D(K+1) = A(K+1) / DEN
          W(K+1) = (Y(K+1)-C(K+1)*W(K+2)) / DEN
   30   Continue
        DEN = B(1) - CRT - C(1) * D(2)
        If (CABS(DEN).NE.0) Then
          Y(1) = (Y(1)-C(1)*W(2)) / DEN
        Else
          Y(1) = (1.,0.)
        End If
        Do 40 J = 2, M
          Y(J) = W(J) - D(J) * Y(J-1)
   40   Continue
      End If
      If (M1.LE.0) Then
        If (M2.LE.0) Go To 60
        RT = BM2(M2)
        M2 = M2 - 1
      Else
        If (M2.LE.0) Then
          RT = BM1(M1)
          M1 = M1 - 1
        Else
          If (ABS(BM1(M1)).GT.ABS(BM2(M2))) Then
            RT = BM1(M1)
            M1 = M1 - 1
          Else
            RT = BM2(M2)
            M2 = M2 - 1
          End If
        End If
      End If
      Y1 = (B(1)-RT) * Y(1) + C(1) * Y(2)
      If (MM.GE.2) Then
c
c matrix multiplication
c
        Do 50 J = 2, MM
          Y2 = A(J) * Y(J-1) + (B(J)-RT) * Y(J) + C(J) * Y(J+1)
          Y(J-1) = Y1
          Y1 = Y2
   50   Continue
      End If
      Y(M) = A(M) * Y(M-1) + (B(M)-RT) * Y(M)
      Y(M-1) = Y1
      IFLG = 1
      Go To 20
   60 If (IA.GT.0) Then
        RT = AA(IA)
        IA = IA - 1
        IFLG = 1
c
c scalar multiplication
c
        Do 70 J = 1, M
          Y(J) = RT * Y(J)
   70   Continue
      End If
      If (IFLG.GT.0) Go To 20
      Do 80 J = 1, M
        YY(J) = REAL(Y(J))
   80 Continue
      Return
      End











      Subroutine CPRODP(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,YY,M,A,B,C,D,U,Y)
c
c prodp applies a sequence of matrix operations to the vector x and
c stores the result in yy       periodic boundary conditions
c and  complex  case
c
c bd,bm1,bm2 are arrays containing roots of certian b polynomials
c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
c aa   array containing scalar multipliers of the vector x
c na is the length of the array aa
c x,yy the matrix operations are applied to x and the result is yy
c a,b,c  are arrays which contain the tridiagonal matrix
c m  is the order of the matrix
c d,u,y are working arrays
c isgn  determines whether or not a change in sign is made
c
      IMPLICIT NONE

      Complex Y, D, U, V, DEN, BH, YM, AM, Y1, Y2, YH, BD, CRT
      Dimension A(1), B(1), C(1), X(1), Y(2), D(1), U(1), BD(1), 
     *    BM1(1), BM2(1), AA(1), YY(1)


      integer J, M, MM, MM2, ND, NM1, NM2, ID, M1, M2, IA
      integer NA, IFLG, K
      real*8 BM1, BM2, AA, X, YY, A, B, C, RT



      Do 10 J = 1, M
        Y(J) = DCMPLX(X(J),0.d0)
   10 Continue
      MM = M - 1
      MM2 = M - 2
      ID = ND
      M1 = NM1
      M2 = NM2
      IA = NA
   20 IFLG = 0
      If (ID.GT.0) Then
        CRT = BD(ID)
        ID = ID - 1
        IFLG = 1
c
c begin solution to system
c
        BH = B(M) - CRT
        YM = Y(M)
        DEN = B(1) - CRT
        D(1) = C(1) / DEN
        U(1) = A(1) / DEN
        Y(1) = Y(1) / DEN
        V = DCMPLX(C(M),0.d0)
        If (MM2.GE.2) Then
          Do 30 J = 2, MM2
            DEN = B(J) - CRT - A(J) * D(J-1)
            D(J) = C(J) / DEN
            U(J) = -A(J) * U(J-1) / DEN
            Y(J) = (Y(J)-A(J)*Y(J-1)) / DEN
            BH = BH - V * U(J-1)
            YM = YM - V * Y(J-1)
            V = -V * D(J-1)
   30     Continue
        End If
        DEN = B(M-1) - CRT - A(M-1) * D(M-2)
        D(M-1) = (C(M-1)-A(M-1)*U(M-2)) / DEN
        Y(M-1) = (Y(M-1)-A(M-1)*Y(M-2)) / DEN
        AM = A(M) - V * D(M-2)
        BH = BH - V * U(M-2)
        YM = YM - V * Y(M-2)
        DEN = BH - AM * D(M-1)
        If (CABS(DEN).NE.0) Then
          Y(M) = (YM-AM*Y(M-1)) / DEN
        Else
          Y(M) = (1.,0.)
        End If
        Y(M-1) = Y(M-1) - D(M-1) * Y(M)
        Do 40 J = 2, MM
          K = M - J
          Y(K) = Y(K) - D(K) * Y(K+1) - U(K) * Y(M)
   40   Continue
      End If
      If (M1.LE.0) Then
        If (M2.LE.0) Go To 60
        RT = BM2(M2)
        M2 = M2 - 1
      Else
        If (M2.LE.0) Then
          RT = BM1(M1)
          M1 = M1 - 1
        Else
          If (DABS(BM1(M1)).GT.DABS(BM2(M2))) Then
            RT = BM1(M1)
            M1 = M1 - 1
          Else
            RT = BM2(M2)
            M2 = M2 - 1
          End If
        End If
      End If
c
c matrix multiplication
c
      YH = Y(1)
      Y1 = (B(1)-RT) * Y(1) + C(1) * Y(2) + A(1) * Y(M)
      If (MM.GE.2) Then
        Do 50 J = 2, MM
          Y2 = A(J) * Y(J-1) + (B(J)-RT) * Y(J) + C(J) * Y(J+1)
          Y(J-1) = Y1
          Y1 = Y2
   50   Continue
      End If
      Y(M) = A(M) * Y(M-1) + (B(M)-RT) * Y(M) + C(M) * YH
      Y(M-1) = Y1
      IFLG = 1
      Go To 20
   60 If (IA.GT.0) Then
        RT = AA(IA)
        IA = IA - 1
        IFLG = 1
c
c scalar multiplication
c
        Do 70 J = 1, M
          Y(J) = RT * Y(J)
   70   Continue
      End If
      If (IFLG.GT.0) Go To 20
      Do 80 J = 1, M
        YY(J) = REAL(Y(J))
   80 Continue
      Return
      End





      Subroutine PPADD(N,IERROR,A,C,CBP,BP,BH)
c
c     ppadd computes the eigenvalues of the periodic tridiagonal matrix
c     with coefficients an,bn,cn
c
c n is the order of the bh and bp polynomials
c     on output bp contians the eigenvalues
c cbp is the same as bp except type complex
c bh is used to temporarily store the roots of the b hat polynomial
c which enters through bp
c
      IMPLICIT NONE
      DOUBLE Complex CX, FSG, HSG, DD
      DOUBLE Complex F, FP, FPP, CDIS, R1, R2, R3, CBP
      Dimension A(1), C(1), BP(1), BH(3), CBP(1)
      real*8 A, C, BP, BH

      real*8 XR, XM



      real*8 SGN

      integer K, NM, NCMPLX, IK
      real*8 EPS, NPP, XL, CNV
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS
      real*8 PSGF, PPSPF, PPSGF

      External PSGF, PPSPF, PPSGF

      real*8 SCNV, DB, BSRH, PSG
      integer IZ, N, IZM, IZM2, J, NT, MODIZ, IS, IF, IG, IT, ICV, I2
      integer I3, NHALF, IERROR


      SCNV = DSQRT(CNV)
      IZ = N
      IZM = IZ - 1
      IZM2 = IZ - 2
      If (BP(N)-BP(1)) 10, 240, 30
   10 Do 20 J = 1, N
        NT = N - J
        BH(J) = BP(NT+1)
   20 Continue
      Go To 50
   30 Do 40 J = 1, N
        BH(J) = BP(J)
   40 Continue
   50 NCMPLX = 0
      MODIZ = MOD(IZ,2)
      IS = 1
      If (MODIZ.NE.0) Then
        If (A(1)) 80, 240, 60
      End If
   60 XL = BH(1)
      DB = BH(3) - BH(1)
   70 XL = XL - DB
      If (PSGF(XL,IZ,C,A,BH).LE.0) Go To 70
      SGN = -1.
      CBP(1) = DCMPLX(BSRH(XL,BH(1),IZ,C,A,BH,PSGF,SGN),0.d0)
      IS = 2
   80 IF = IZ - 1
      If (MODIZ.NE.0) Then
        If (A(1)) 90, 240, 110
      End If
   90 XR = BH(IZ)
      DB = BH(IZ) - BH(IZ-2)
  100 XR = XR + DB
      If (PSGF(XR,IZ,C,A,BH).LT.0) Go To 100
      SGN = 1.
      CBP(IZ) = DCMPLX(BSRH(BH(IZ),XR,IZ,C,A,BH,PSGF,SGN),0.d0)
      IF = IZ - 2
  110 Do 180 IG = IS, IF, 2
        XL = BH(IG)
        XR = BH(IG+1)
        SGN = -1.
        XM = BSRH(XL,XR,IZ,C,A,BH,PPSPF,SGN)
        PSG = PSGF(XM,IZ,C,A,BH)
        If (ABS(PSG).GT.EPS) Then
          If (PSG*PPSGF(XM,IZ,C,A,BH)) 120, 130, 140
c
c     case of a real zero
c
  120     SGN = 1.
          CBP(IG) = DCMPLX(BSRH(BH(IG),XM,IZ,C,A,BH,PSGF,SGN),0.d0)
          SGN = -1.
          CBP(IG+1) = DCMPLX(BSRH(XM,BH(IG+1),IZ,C,A,BH,PSGF,SGN),0.d0)
          Go To 180
        End If
c
c     case of a multiple zero
c
  130   CBP(IG) = DCMPLX(XM,0.d0)
        CBP(IG+1) = DCMPLX(XM,0.d0)
        Go To 180
c
c     case of a complex zero
c
  140   IT = 0
        ICV = 0
        CX = DCMPLX(XM,0.d0)
  150   FSG = (1.,0.)
        HSG = (1.,0.)
        FP = (0.,0.)
        FPP = (0.,0.)
        Do 160 J = 1, IZ
          DD = 1. / (CX-BH(J))
          FSG = FSG * A(J) * DD
          HSG = HSG * C(J) * DD
          FP = FP + DD
          FPP = FPP - DD * DD
  160   Continue
        If (MODIZ.EQ.0) Then
          F = (1.,0.) - FSG - HSG
        Else
          F = (1.,0.) + FSG + HSG
        End If
        I3 = 0
        If (CDABS(FP).GT.0) Then
          I3 = 1
          R3 = -F / FP
        End If
        I2 = 0
        If (CDABS(FPP).GT.0) Then
          I2 = 1
          CDIS = CDSQRT(FP**2-2.*F*FPP)
          R1 = CDIS - FP
          R2 = -FP - CDIS
          If (CDABS(R1).GT.CDABS(R2)) Then
            R1 = R1 / FPP
          Else
            R1 = R2 / FPP
          End If
          R2 = 2. * F / FPP / R1
          If (CDABS(R2).LT.CDABS(R1)) R1 = R2
          If (I3.LE.0) Go To 170
          If (CDABS(R3).LT.CDABS(R1)) R1 = R3
        Else
          R1 = R3
        End If
  170   CX = CX + R1
        IT = IT + 1
        If (IT.GT.50) Go To 240
        If (CDABS(R1).GT.SCNV) Go To 150
        If (ICV.LE.0) Then
          ICV = 1
          Go To 150
        End If
        CBP(IG) = CX
        CBP(IG+1) = CONJG(CX)
  180 Continue
      If (CDABS(CBP(N))-CDABS(CBP(1))) 190, 240, 210
  190 NHALF = N / 2
      Do 200 J = 1, NHALF
        NT = N - J
        CX = CBP(J)
        CBP(J) = CBP(NT+1)
        CBP(NT+1) = CX
  200 Continue
  210 NCMPLX = 1
      Do 220 J = 2, IZ
        If (DIMAG(CBP(J)).NE.0) Go To 250
  220 Continue
      NCMPLX = 0
      Do 230 J = 2, IZ
        BP(J) = DREAL(CBP(J))
  230 Continue
      Go To 250
  240 IERROR = 4
  250 Continue
      Return
      End





      Function PSGF(X,IZ,C,A,BH)
      IMPLICIT NONE
      DOUBLE PRECISION PSGF
      Dimension A(1), C(1), BH(1)
      real*8 A, X, C, BH
      integer IZ

      real*8 FSG, HSG, DD
      integer J

      FSG = 1.
      HSG = 1.
      Do 10 J = 1, IZ
        DD = 1. / (X-BH(J))
        FSG = FSG * A(J) * DD
        HSG = HSG * C(J) * DD
   10 Continue
      If (MOD(IZ,2).EQ.0) Then
        PSGF = 1. - FSG - HSG
        Return
      End If
      PSGF = 1. + FSG + HSG
      Return
      End










      Function BSRH(XLL,XRR,IZ,C,A,BH,F,SGN)
      IMPLICIT NONE
      DOUBLE PRECISION BSRH
      DOUBLE PRECISION F
      Dimension A(1), C(1), BH(1)
      real*8 XLL, XRR, A, C, BH, SGN
      integer IZ

      integer K, NM, NCMPLX, IK
      real*8 EPS, NPP, XL, CNV
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS

      real*8 DX, X, XR

      XL = XLL
      XR = XRR
      DX = .5 * ABS(XR-XL)
   10 X = .5 * (XL+XR)
      If (SGN*F(X,IZ,C,A,BH)) 30, 50, 20
   20 XR = X
      Go To 40
   30 XL = X
   40 DX = .5 * DX
      If (DX.GT.CNV) Go To 10
   50 BSRH = .5 * (XL+XR)
      Return
      End








      Function PPSGF(X,IZ,C,A,BH)
      IMPLICIT NONE
      real*8 PPSGF
      Dimension A(1), C(1), BH(1)
      real*8 A, C, BH, X
      integer J, IZ

      real*8 SUM
      SUM = 0.
      Do 10 J = 1, IZ
        SUM = SUM - 1. / (X-BH(J)) ** 2
   10 Continue
      PPSGF = SUM
      Return
      End







      Function PPSPF(X,IZ,C,A,BH)
      IMPLICIT NONE
      real*8 PPSPF
      Dimension A(1), C(1), BH(1)
      real*8 SUM, A, C, BH, X
      integer J, IZ
      SUM = 0.
      Do 10 J = 1, IZ
        SUM = SUM + 1. / (X-BH(J))
   10 Continue
      PPSPF = SUM
      Return
      End








      Subroutine COMPB(N,IERROR,AN,BN,CN,B,AH,BH)
c
c     compb computes the roots of the b polynomials using subroutine
c     tevls which is a modification the eispack program tqlrat.
c     ierror is set to 4 if either tevls fails or if a(j+1)*c(j) is
c     less than zero for some j.  ah,bh are temporary work arrays.
c
      IMPLICIT NONE
      Dimension AN(1), BN(1), CN(1), B(1), AH(1), BH(1)
      real*8 AN, BN, CN, B, AH, BH
      integer N, IERROR
      real*8 BNORM, EPMACH


 

      integer J, K, NM, NCMPLX, IK, DUM, IF, KDO, L, IR, I2, I4, IPL 
      integer IFD, I, IB, NB, JS, JF, LS, LH, NMP, L1, L2, J2, J1, N2M2
      real*8 EPS, NPP, CNV, ARG, D1, D2, D3
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, EPS

      EPS = EPMACH(DUM)
      BNORM = DABS(BN(1))
      Do 10 J = 2, NM
        BNORM = DMAX1(BNORM,DABS(BN(J)))
        ARG = AN(J) * CN(J-1)
        If (ARG.LT.0) Go To 110
        B(J) = SIGN(SQRT(ARG),AN(J))
   10 Continue
      CNV = EPS * BNORM
      IF = 2 ** K
      KDO = K - 1
      Do 50 L = 1, KDO
        IR = L - 1
        I2 = 2 ** IR
        I4 = I2 + I2
        IPL = I4 - 1
        IFD = IF - I4
        Do 40 I = I4, IFD, I4
          Call INDXB(I,L,IB,NB)
          If (NB.LE.0) Go To 50
          JS = I - IPL
          JF = JS + NB - 1
          LS = 0
          Do 20 J = JS, JF
            LS = LS + 1
            BH(LS) = BN(J)
            AH(LS) = B(J)
   20     Continue
          Call TEVLS(NB,BH,AH,IERROR)
          If (IERROR.NE.0) Go To 100
          LH = IB - 1
          Do 30 J = 1, NB
            LH = LH + 1
            B(LH) = -BH(J)
   30     Continue
   40   Continue
   50 Continue
      Do 60 J = 1, NM
        B(J) = -BN(J)
   60 Continue
      If (NPP.EQ.0) Then
        NMP = NM + 1
        NB = NM + NMP
        Do 70 J = 1, NB
          L1 = MOD(J-1,NMP) + 1
          L2 = MOD(J+NM-1,NMP) + 1
          ARG = AN(L1) * CN(L2)
          If (ARG.LT.0) Go To 110
          BH(J) = SIGN(SQRT(ARG),-AN(L1))
          AH(J) = -BN(L1)
   70   Continue
        Call TEVLS(NB,AH,BH,IERROR)
        If (IERROR.NE.0) Go To 100
        Call INDXB(IF,K-1,J2,LH)
        Call INDXB(IF/2,K-1,J1,LH)
        J2 = J2 + 1
        LH = J2
        N2M2 = J2 + NM + NM - 2
   80   D1 = DABS(B(J1)-B(J2-1))
        D2 = DABS(B(J1)-B(J2))
        D3 = DABS(B(J1)-B(J2+1))
        If (.NOT.((D2.LT.D1).AND.(D2.LT.D3))) Then
          B(LH) = B(J2)
          J2 = J2 + 1
          LH = LH + 1
          If (J2-N2M2) 80, 80, 90
        End If
        J2 = J2 + 1
        J1 = J1 + 1
        If (J2.LE.N2M2) Go To 80
   90   B(LH) = B(N2M2+1)
        Call INDXB(IF,K-1,J1,J2)
        J2 = J1 + NMP + NMP
        Call PPADD(NM+1,IERROR,AN,CN,B(J1),B(J1),B(J2))
      End If
      Return
  100 IERROR = 4
      Return
  110 IERROR = 5
      Return
      End








      Subroutine TEVLS(N,D,E2,IERR)
c
      IMPLICIT NONE
      Integer I, J, L, M, N, II, L1, MML, IERR
      Real*8 D(N), E2(N)
      Real*8 B, C, F, G, H, P, R, S
c
c     real sqrt,abs,sign
c

      integer K, NM, NCMPLX, IK, NHALF, NTOP
      real*8 MACHEP, NPP, CNV, DHOLD
      Common /CBLKT/ K, NM, NCMPLX, IK, NPP, CNV, MACHEP

c
c     this subroutine is a modification of the eispack subroutine tqlrat
c     algorithm 464, comm. acm 16, 689(1973) by reinsch.
c
c     this subroutine finds the eigenvalues of a symmetric
c     tridiagonal matrix by the rational ql method.
c
c     on input-
c
c        n is the order of the matrix,
c
c        d contains the diagonal elements of the input matrix,
c
c        e2 contains the                subdiagonal elements of the
c          input matrix in its last n-1 positions.  e2(1) is arbitrary.
c
c      on output-
c
c        d contains the eigenvalues in ascending order.  if an
c          error exit is made, the eigenvalues are correct and
c          ordered for indices 1,2,...ierr-1, but may not be
c          the smallest eigenvalues,
c
c        e2 has been destroyed,
c
c        ierr is set to
c          zero       for normal return,
c          j          if the j-th eigenvalue has not been
c                     determined after 30 iterations.
c
c     questions and comments should be directed to b. s. garbow,
c     applied mathematics division, argonne national laboratory
c
c
c     ********** machep is a machine dependent parameter specifying
c                the relative precision of floating point arithmetic.
c
c                **********
c
      IERR = 0
      If (N.NE.1) Then
c
        Do 10 I = 2, N
          E2(I-1) = E2(I) * E2(I)
   10   Continue
c
        F = 0.0
        B = 0.0
        E2(N) = 0.0
c
        Do 90 L = 1, N
          J = 0
          H = MACHEP * (DABS(D(L))+DSQRT(E2(L)))
          If (B.LE.H) Then
            B = H
            C = B * B
          End If
c
c     ********** look for small squared sub-diagonal element **********
c
          Do 20 M = L, N
            If (E2(M).LE.C) Go To 30
c
c     ********** e2(n) is always zero, so there is no exit
c                through the bottom of the loop **********
c
   20     Continue
c
   30     If (M.NE.L) Then
   40       If (J.EQ.30) Go To 110
            J = J + 1
c
c     ********** form shift **********
c
            L1 = L + 1
            S = SQRT(E2(L))
            G = D(L)
            P = (D(L1)-G) / (2.0*S)
            R = SQRT(P*P+1.0)
            D(L) = S / (P+SIGN(R,P))
            H = G - D(L)
c
            Do 50 I = L1, N
              D(I) = D(I) - H
   50       Continue
c
            F = F + H
c
c     ********** rational ql transformation **********
c
            G = D(M)
            If (G.EQ.0.0) G = B
            H = G
            S = 0.0
            MML = M - L
c
c     ********** for i=m-1 step -1 until l do -- **********
c
            Do 60 II = 1, MML
              I = M - II
              P = G * H
              R = P + E2(I)
              E2(I+1) = S * R
              S = E2(I) / R
              D(I+1) = H + S * (H+D(I))
              G = D(I) - E2(I) / G
              If (G.EQ.0.0) G = B
              H = G * P / R
   60       Continue
c
            E2(L) = S * G
            D(L) = H
c
c     ********** guard against underflowed h **********
c
            If (H.NE.0.0) Then
              If (ABS(E2(L)).GT.ABS(C/H)) Then
                E2(L) = H * E2(L)
                If (E2(L).NE.0.0) Go To 40
              End If
            End If
          End If
          P = D(L) + F
c
c     ********** order eigenvalues **********
c
          If (L.NE.1) Then
c
c     ********** for i=l step -1 until 2 do -- **********
c
            Do 70 II = 2, L
              I = L + 2 - II
              If (P.GE.D(I-1)) Go To 80
              D(I) = D(I-1)
   70       Continue
          End If
c
          I = 1
   80     D(I) = P
   90   Continue
c
        If (ABS(D(N)).GE.ABS(D(1))) Go To 120
        NHALF = N / 2
        Do 100 I = 1, NHALF
          NTOP = N - I
          DHOLD = D(I)
          D(I) = D(NTOP+1)
          D(NTOP+1) = DHOLD
  100   Continue
        Go To 120
c
c     ********** set error -- no convergence to an
c                eigenvalue after 30 iterations **********
c
  110   IERR = L
      End If
  120 Return
c
c     ********** last card of tqlrat **********
c
c
c revision history---
c
c december 1979    first added to nssl
c-----------------------------------------------------------------------
      End
c package comf           the entries in this package are lowlevel
c                        entries, supporting ulib packages blktri
c                        and cblktri. that is, these routines are
c                        not called directly by users, but rather
c                        by entries within blktri and cblktri.
c                        description of entries epmach and pimach
c                        follow below.
c
c latest revision        january 1985
c
c special conditions     none
c
c i/o                    none
c
c precision              single
c
c required library       none
c files
c
c language               fortran
c ********************************************************************
c
c function epmach (dum)
c
c purpose                to compute an approximate machine accuracy
c                        epsilon according to the following definition:
c                        epsilon is the smallest number such that
c                        (1.+epsilon).gt.1.)
c
c usage                  eps = epmach (dum)
c
c arguments
c on input               dum
c                          dummy value
c
c arguments
c on output              none
c
c history                the original version, written when the
c                        blktri package was converted from the
c                        cdc 7600 to run on the cray-1, calculated
c                        machine accuracy by successive divisions
c                        by 10.  use of this constant caused blktri
c                        to compute solutions on the cray-1 with four
c                        fewer places of accuracy than the version
c                        on the 7600.  it was found that computing
c                        machine accuracy by successive divisions
c                        of 2 produced a machine accuracy 29% less
c                        than the value generated by successive
c                        divisions by 10, and that use of this
c                        machine constant in the blktri package
c                        recovered the accuracy that appeared to
c                        be lost on conversion.
c
c algorithm              computes machine accuracy by successive
c                        divisions of two.
c
c portability            this code will execute on machines other
c                        than the cray1, but the returned value may
c                        be unsatisfactory.  see history above.
c ********************************************************************
c
c function pimach (dum)
c
c purpose                to supply the value of the constant pi
c                        correct to machine precision where
c                        pi=3.141592653589793238462643383279502884197
c                             1693993751058209749446
c
c usage                  pi = pimach (dum)
c
c arguments
c on input               dum
c                          dummy value
c
c arguments
c on output              none
c
c algorithm              the value of pi is set in a constant.
c
c portability            this entry is portable, but users should
c                        check to see whether greater accuracy is
c                        required.
c
c***********************************************************************
      Function EPMACH(DUM)
      IMPLICIT NONE
      real*8 EPMACH, V, EPS
      integer DUM
      Common /VALUE/ V
      EPS = 1.
   10 EPS = EPS / 2.
      Call STORE(EPS+1.)
      If (V.GT.1.) Go To 10
      EPMACH = 100. * EPS
      Return
      End
      Subroutine STORE(X)
      IMPLICIT NONE
      real*8 X, V
      Common /VALUE/ V
      V = X
      Return
      End







C      Function PIMACH(DUM)
C      IMPLICIT NONE
C      real*8 PIMACH
c     pi=3.1415926535897932384626433832795028841971693993751058209749446
c
C      PIMACH = 3.14159265358979
C      Return
C      End
c package sepaux         contains no user entry points.
c
c latest revision        march 1985
c
c purpose                this package contains auxiliary routines for
c                        ncar public software packages such as sepeli
c                        and sepx4.
c
c usage                  since this package contains no user entries,
c                        no usage instructions or argument descriptions
c                        are given here.
c
c special conditions     none
c
c i/o                    none
c
c precision              single
c
c required library       none
c files
c
c language               fortran
c
c history                developed in the late 1970's by john c. adams
c                        of ncar's scienttific computing division.
c
c portability            fortran 66
c **********************************************************************
      Subroutine SEPORT(USOL,IDMN,ZN,ZM,PERTRB)
c
c     this subroutine orthoganalizes the array usol with respect to
c     the constant array in a weighted least squares norm
c
      IMPLICIT NONE
      integer KSWX, KSWY, K , L
      integer AIT, BIT, CIT, DIT, IS, MIT, NIT
      real*8 MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4
      Common /SPLP/ KSWX, KSWY, K, L, 
     *    MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4,
     *    AIT, BIT, CIT, DIT, MIT, NIT, IS
      integer IDMN
      Dimension USOL(IDMN,1), ZN(1), ZM(1)
      real*8 USOL, ZN, ZM, UTE, ETE, PERTRB
      integer ISTR, JSTR, IFNL, JFNL, I, II, J, JJ


      ISTR = IS
      IFNL = MS
      JSTR = JS
      JFNL = NS
c
c     compute weighted inner products
c
      UTE = 0.0
      ETE = 0.0
      Do 20 I = IS, MS
        II = I - IS + 1
        Do 10 J = JS, NS
          JJ = J - JS + 1
          ETE = ETE + ZM(II) * ZN(JJ)
          UTE = UTE + USOL(I,J) * ZM(II) * ZN(JJ)
   10   Continue
   20 Continue
c
c     set perturbation parameter
c
      PERTRB = UTE / ETE
c
c     subtract off constant pertrb
c
      Do 40 I = ISTR, IFNL
        Do 30 J = JSTR, JFNL
          USOL(I,J) = USOL(I,J) - PERTRB
   30   Continue
   40 Continue
      Return
      End








      Subroutine SEPMIN(USOL,IDMN,ZN,ZM,PERTB)
c
c     this subroutine orhtogonalizes the array usol with respect to
c     the constant array in a weighted least squares norm
c
      IMPLICIT NONE
      integer KSWX, KSWY, K, IDMN, ISTR, JSTR,L, IFNL, JFNL
      integer AIT, BIT, CIT, DIT, IS, MIT, NIT
      real*8 MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4
      real*8 UTE, ETE, PERTB, PERTRB
      integer I, II, J, JJ
      Common /SPLP/ KSWX, KSWY, K, L, 
     *    MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4,
     *    AIT, BIT, CIT, DIT, MIT, NIT, IS
      Dimension USOL(IDMN,1), ZN(1), ZM(1)
      real*8 USOL, ZN, ZM
c
c     entry at sepmin occurrs when the final solution is
c     to be minimized with respect to the weighted
c     least squares norm
c
      ISTR = 1
      IFNL = K
      JSTR = 1
      JFNL = L
c
c     compute weighted inner products
c
      UTE = 0.0
      ETE = 0.0
      Do 20 I = IS, MS
        II = I - IS + 1
        Do 10 J = JS, NS
          JJ = J - JS + 1
          ETE = ETE + ZM(II) * ZN(JJ)
          UTE = UTE + USOL(I,J) * ZM(II) * ZN(JJ)
   10   Continue
   20 Continue
c
c     set perturbation parameter
c
      PERTRB = UTE / ETE
c
c     subtract off constant pertrb
c
      Do 40 I = ISTR, IFNL
        Do 30 J = JSTR, JFNL
          USOL(I,J) = USOL(I,J) - PERTRB
   30   Continue
   40 Continue
      Return
      End









      Subroutine SEPTRI(N,A,B,C,D,U,Z)
c
c     this subroutine solves for a non-zero eigenvector corresponding
c     to the zero eigenvalue of the transpose of the rank
c     deficient one matrix with subdiagonal a, diagonal b, and
c     superdiagonal c , with a(1) in the (1,n) position, with
c     c(n) in the (n,1) position, and all other elements zero.
c
      IMPLICIT NONE
      integer N
      Dimension A(N), B(N), C(N), D(N), U(N), Z(N)
      real*8 A, B, C, D, U, Z

      real*8 BN, V, DEN, AN
      integer K, NM1, NM2, J
      BN = B(N)
      D(1) = A(2) / B(1)
      V = A(1)
      U(1) = C(N) / B(1)
      NM2 = N - 2
      Do 10 J = 2, NM2
        DEN = B(J) - C(J-1) * D(J-1)
        D(J) = A(J+1) / DEN
        U(J) = -C(J-1) * U(J-1) / DEN
        BN = BN - V * U(J-1)
        V = -V * D(J-1)
   10 Continue
      DEN = B(N-1) - C(N-2) * D(N-2)
      D(N-1) = (A(N)-C(N-2)*U(N-2)) / DEN
      AN = C(N-1) - V * D(N-2)
      BN = BN - V * U(N-2)
      DEN = BN - AN * D(N-1)
c
c     set last component equal to one
c
      Z(N) = 1.0
      Z(N-1) = -D(N-1)
      NM1 = N - 1
      Do 20 J = 2, NM1
        K = N - J
        Z(K) = -D(K) * Z(K+1) - U(K) * Z(N)
   20 Continue
      Return
      End

























      Subroutine SEPDX(U,IDMN,I,J,UXXX,UXXXX)
c
c     this program computes second order finite difference
c     approximations to the third and fourth x
c     partial derivatives of u at the (i,j) mesh point
c
      IMPLICIT NONE
      integer KSWX, KSWY, K, L, I, J
      integer AIT, BIT, CIT, DIT, IS, MIT, NIT
      real*8 MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4
      Common /SPLP/ KSWX, KSWY, K, L, 
     *    MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4,
     *    AIT, BIT, CIT, DIT, MIT, NIT, IS
      integer IDMN
      Dimension U(IDMN,1)
      real*8 U, UXXX, UXXXX
      If (I.LE.2.OR.I.GE.(K-1)) Then
        If (I.NE.1) Then
          If (I.EQ.2) Go To 10
          If (I.EQ.K-1) Go To 20
          If (I.EQ.K) Go To 30
        End If
c
c     compute partial derivative approximations at x=a
c
        If (KSWX.NE.1) Then
          UXXX = (-5.0*U(1,J)+18.0*U(2,J)-24.0*U(3,J)+14.0*U(4,J)-3.0*
     *        U(5,J)) / (TDLX3)
          UXXXX = (3.0*U(1,J)-14.0*U(2,J)+26.0*U(3,J)-24.0*U(4,J)+11.0*
     *        U(5,J)-2.0*U(6,J)) / DLX4
          Return
        End If
c
c     periodic at x=a
c
        UXXX = (-U(K-2,J)+2.0*U(K-1,J)-2.0*U(2,J)+U(3,J)) / (TDLX3)
        UXXXX = (U(K-2,J)-4.0*U(K-1,J)+6.0*U(1,J)-4.0*U(2,J)+U(3,J)) / 
     *      DLX4
        Return
c
c     compute partial derivative approximations at x=a+dlx
c
   10   If (KSWX.NE.1) Then
          UXXX = (-3.0*U(1,J)+10.0*U(2,J)-12.0*U(3,J)+6.0*U(4,J)-
     *        U(5,J)) / TDLX3
          UXXXX = (2.0*U(1,J)-9.0*U(2,J)+16.0*U(3,J)-14.0*U(4,J)+6.0*
     *        U(5,J)-U(6,J)) / DLX4
          Return
        End If
c
c     periodic at x=a+dlx
c
        UXXX = (-U(K-1,J)+2.0*U(1,J)-2.0*U(3,J)+U(4,J)) / (TDLX3)
        UXXXX = (U(K-1,J)-4.0*U(1,J)+6.0*U(2,J)-4.0*U(3,J)+U(4,J)) / 
     *      DLX4
        Return
      End If
c
c     compute partial derivative approximations on the interior
c
      UXXX = (-U(I-2,J)+2.0*U(I-1,J)-2.0*U(I+1,J)+U(I+2,J)) / TDLX3
      UXXXX = (U(I-2,J)-4.0*U(I-1,J)+6.0*U(I,J)-4.0*U(I+1,J)+U(I+2,J)) /
     *    DLX4
      Return
c
c     compute partial derivative approximations at x=b-dlx
c
   20 If (KSWX.NE.1) Then
        UXXX = (U(K-4,J)-6.0*U(K-3,J)+12.0*U(K-2,J)-10.0*U(K-1,J)+3.0*
     *      U(K,J)) / TDLX3
        UXXXX = (-U(K-5,J)+6.0*U(K-4,J)-14.0*U(K-3,J)+16.0*U(K-2,J)-9.0*
     *      U(K-1,J)+2.0*U(K,J)) / DLX4
        Return
      End If
c
c     periodic at x=b-dlx
c
      UXXX = (-U(K-3,J)+2.0*U(K-2,J)-2.0*U(1,J)+U(2,J)) / TDLX3
      UXXXX = (U(K-3,J)-4.0*U(K-2,J)+6.0*U(K-1,J)-4.0*U(1,J)+U(2,J)) / 
     *    DLX4
      Return
c
c     compute partial derivative approximations at x=b
c
   30 UXXX = -(3.0*U(K-4,J)-14.0*U(K-3,J)+24.0*U(K-2,J)-18.0*U(K-1,J)+
     *    5.0*U(K,J)) / TDLX3
      UXXXX = (-2.0*U(K-5,J)+11.0*U(K-4,J)-24.0*U(K-3,J)+26.0*U(K-2,J)-
     *    14.0*U(K-1,J)+3.0*U(K,J)) / DLX4
      Return
      End






















      Subroutine SEPDY(U,IDMN,I,J,UYYY,UYYYY)
      IMPLICIT NONE
c
c     this program computes second order finite difference
c     approximations to the third and fourth y
c     partial derivatives of u at the (i,j) mesh point
c
      integer KSWX, KSWY, K 
      integer AIT, BIT, CIT, DIT, IS, MIT, NIT
      real*8 MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4
      Common /SPLP/ KSWX, KSWY, K, L, 
     *    MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4,
     *    AIT, BIT, CIT, DIT, MIT, NIT, IS
      integer IDMN, I, J, L
      Dimension U(IDMN,IDMN)
      real*8 U, UYYY, UYYYY

      If (J.LE.2.OR.J.GE.(L-1)) Then
        If (J.NE.1) Then
          If (J.EQ.2) Go To 10
          If (J.EQ.L-1) Go To 20
          If (J.EQ.L) Go To 30
        End If
c
c     compute partial derivative approximations at y=c
c
        If (KSWY.NE.1) Then
          UYYY = (-5.0*U(I,1)+18.0*U(I,2)-24.0*U(I,3)+14.0*U(I,4)-3.0*
     *        U(I,5)) / TDLY3
          UYYYY = (3.0*U(I,1)-14.0*U(I,2)+26.0*U(I,3)-24.0*U(I,4)+11.0*
     *        U(I,5)-2.0*U(I,6)) / DLY4
          Return
        End If
c
c     periodic at x=a
c
        UYYY = (-U(I,L-2)+2.0*U(I,L-1)-2.0*U(I,2)+U(I,3)) / TDLY3
        UYYYY = (U(I,L-2)-4.0*U(I,L-1)+6.0*U(I,1)-4.0*U(I,2)+U(I,3)) / 
     *      DLY4
        Return
c
c     compute partial derivative approximations at y=c+dly
c
   10   If (KSWY.NE.1) Then
          UYYY = (-3.0*U(I,1)+10.0*U(I,2)-12.0*U(I,3)+6.0*U(I,4)-
     *        U(I,5)) / TDLY3
          UYYYY = (2.0*U(I,1)-9.0*U(I,2)+16.0*U(I,3)-14.0*U(I,4)+6.0*
     *        U(I,5)-U(I,6)) / DLY4
          Return
        End If
c
c     periodic at y=c+dly
c
        UYYY = (-U(I,L-1)+2.0*U(I,1)-2.0*U(I,3)+U(I,4)) / TDLY3
        UYYYY = (U(I,L-1)-4.0*U(I,1)+6.0*U(I,2)-4.0*U(I,3)+U(I,4)) / 
     *      DLY4
        Return
      End If
c
c     compute partial derivative approximations on the interior
c
      UYYY = (-U(I,J-2)+2.0*U(I,J-1)-2.0*U(I,J+1)+U(I,J+2)) / TDLY3
      UYYYY = (U(I,J-2)-4.0*U(I,J-1)+6.0*U(I,J)-4.0*U(I,J+1)+U(I,J+2)) /
     *    DLY4
      Return
c
c     compute partial derivative approximations at y=d-dly
c
   20 If (KSWY.NE.1) Then
        UYYY = (U(I,L-4)-6.0*U(I,L-3)+12.0*U(I,L-2)-10.0*U(I,L-1)+3.0*
     *      U(I,L)) / TDLY3
        UYYYY = (-U(I,L-5)+6.0*U(I,L-4)-14.0*U(I,L-3)+16.0*U(I,L-2)-9.0*
     *      U(I,L-1)+2.0*U(I,L)) / DLY4
        Return
      End If
c
c     periodic at y=d-dly
c
      UYYY = (-U(I,L-3)+2.0*U(I,L-2)-2.0*U(I,1)+U(I,2)) / TDLY3
      UYYYY = (U(I,L-3)-4.0*U(I,L-2)+6.0*U(I,L-1)-4.0*U(I,1)+U(I,2)) / 
     *    DLY4
      Return
c
c     compute partial derivative approximations at y=d
c
   30 UYYY = -(3.0*U(I,L-4)-14.0*U(I,L-3)+24.0*U(I,L-2)-18.0*U(I,L-1)+
     *    5.0*U(I,L)) / TDLY3
      UYYYY = (-2.0*U(I,L-5)+11.0*U(I,L-4)-24.0*U(I,L-3)+26.0*U(I,L-2)-
     *    14.0*U(I,L-1)+3.0*U(I,L)) / DLY4
      Return
c
c revision history---
c
c december 1979    first added to nssl
c-----------------------------------------------------------------------
      End






      Subroutine COFX(XX,AFUN,BFUN,CFUN)
c
c     subroutine to compute the coefficients of the elliptic equation
c     solved to fill in the grid.  it is passed (along with cofy) to
c     the elliptic solver sepeli.
      Include 'mexsepeli.inc'
      integer I
      real*8 DXDI, AFUN, BFUN, CFUN, XX
      I = DNINT(XX)
      DXDI = 0.5 * (SXI(I+1)-SXI(I-1))
      AFUN = 1. / DXDI ** 2
      BFUN = (-SXI(I+1)+2.*SXI(I)-SXI(I-1)) / DXDI ** 3
      CFUN = 0.
      Return
      End




      Subroutine COFY(YY,DFUN,EFUN,FFUN)
      Include 'mexsepeli.inc'
      integer J
      real*8 DEDJ, DFUN, EFUN, FFUN, YY
      J = DNINT(YY)
      DEDJ = 0.5 * (SETA(J+1)-SETA(J-1))
      DFUN = 1. / DEDJ ** 2
      EFUN = (-SETA(J+1)+2.*SETA(J)-SETA(J-1)) / DEDJ ** 3
      FFUN = 0.
      Return
      End




