C      ALGORITHM 732, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 20, NO. 3, SEPTEMBER, 1994, PP. 247-261.
cat > Readme <<C-END-OF-FILE
       This file contains instructions on how to produce and run
       the components of a FORTRAN software package of routines to
       solve finite difference approximations of certain self-adjoint
       elliptic equations in two-dimensional irregular domains.

	Descriptions of the methodology may be found in the article
	'Solvers for self-adjoint elliptic problems in irregular
 	two-dimensional domains' by P.Cummins and G.Vallis
	and in 'A fast semi--direct method for the numerical
	solution of non--separable elliptic equations in
	irregular domains. J. Comp. Phys. 82, 1989, 398--412.
	by A Pares-Sierra and G. Vallis, and references therein.

       The package is distributed as a single file called Solvers.
       To split Solvers into its components the following UNIX
       command should be executed:

                sh < Solvers

       This should produce the following files:

                Readme
                Makefile
                solvers.f
                poiss.f
                linpack.f
                testcapc.f
                testdcapc.f
                testcapcn.f
                testdcapcn.f
                testreccn.f
                testdreccn.f

       To compile and link files and then run executable code, the
       following UNIX command should be issued:

                       make

       The Makefile assumes the user's FORTRAN compiler is invoked
       with f77 and accepts a -O optimizing directive. If this
       is not the case, then the FC and FFLAGS macros in the
       Makefile should be redefined appropriately.

       The output from the execution of the test programs is found
       in the *.out files which are created at execution time.


       To produce executable code for an individual test program,
       a UNIX command such as the following should be issued:

                 make testcapc

       To clean up file space and to delete object and executable
       files, the following UNIX command should be issued:

                 make clean
C-END-OF-FILE
cat > Makefile <<C-END-OF-FILE
#  Makefile for elliptic solver routines

FC = f77
FFLAGS = -O

ALL = testcapc testdcapc testcapcn testdcapcn testreccn testdreccn

runall : $(ALL)
	testcapc > testcapc.out
	testdcapc > testdcapc.out
	testcapcn > testcapcn.out
	testdcapcn > testdcapcn.out
	testreccn > testreccn.out
	testdreccn > testdreccn.out

#  testcapc

OBJECTS1= testcapc.o solvers.o poiss.o linpack.o
testcapc: $(OBJECTS1)
	$(FC) -o testcapc $(OBJECTS1)

#  testdcapc

OBJECTS2= testdcapc.o solvers.o poiss.o linpack.o
testdcapc: $(OBJECTS2)
	$(FC) -o testdcapc $(OBJECTS2)

#  testcapcn

OBJECTS3= testcapcn.o solvers.o poiss.o linpack.o
testcapcn: $(OBJECTS3)
	$(FC) -o testcapcn $(OBJECTS3)

#  testdcapcn

OBJECTS4= testdcapcn.o solvers.o poiss.o linpack.o
testdcapcn: $(OBJECTS4)
	$(FC) -o testdcapcn $(OBJECTS4)

#  testreccn

OBJECTS5= testreccn.o solvers.o poiss.o linpack.o
testreccn: $(OBJECTS5)
	$(FC) -o testreccn $(OBJECTS5)

#  testdreccn

OBJECTS6= testdreccn.o solvers.o poiss.o linpack.o
testdreccn: $(OBJECTS6)
	$(FC) -o testdreccn $(OBJECTS6)

clean :
	rm -f core $(ALL) *.o
C-END-OF-FILE
cat > linpack.f <<C-END-OF-FILE
      subroutine dgefa(a,lda,n,ipvt,info)
      integer lda,n,ipvt(1),info
      double precision a(lda,1)
c
c     dgefa factors a double precision matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c
c     internal variables
c
      double precision t
      integer idamax,j,k,kp1,l,nm1
c
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0d0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -1.0d0/a(k,k)
            call dscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0d0) info = n
      return
      end
      subroutine daxpy(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
      subroutine  dscal(n,da,dx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified to correct problem with negative increment, 8/21/90.
c
      double precision da,dx(1)
      integer i,incx,ix,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      do 10 i = 1,n
        dx(ix) = da*dx(ix)
        ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end
      integer function idamax(n,dx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c     modified to correct problem with negative increment, 8/21/90.
c
      double precision dx(1),dmax
      integer i,incx,ix,n
c
      idamax = 0
      if( n .lt. 1 ) return
      idamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      dmax = dabs(dx(ix))
      ix = ix + incx
      do 10 i = 2,n
         if(dabs(dx(ix)).le.dmax) go to 5
         idamax = i
         dmax = dabs(dx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 dmax = dabs(dx(1))
      do 30 i = 2,n
         if(dabs(dx(i)).le.dmax) go to 30
         idamax = i
         dmax = dabs(dx(i))
   30 continue
      return
      end
      subroutine sgefa(a,lda,n,ipvt,info)
      integer lda,n,ipvt(1),info
      real a(lda,1)
c
c     sgefa factors a real matrix by gaussian elimination.
c
c     sgefa is usually called by sgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for sgeco) = (1 + 9/n)*(time for sgefa) .
c
c     on entry
c
c        a       real(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that sgesl or sgedi will divide by zero
c                     if called.  use  rcond  in sgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas saxpy,sscal,isamax
c
c     internal variables
c
      real t
      integer isamax,j,k,kp1,l,nm1
c
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = isamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0e0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -1.0e0/a(k,k)
            call sscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0e0) info = n
      return
      end
      subroutine saxpy(n,sa,sx,incx,sy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),sa
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        sy(iy) = sy(iy) + sa*sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sy(i) = sy(i) + sa*sx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        sy(i) = sy(i) + sa*sx(i)
        sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
        sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
        sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
   50 continue
      return
      end
      subroutine sscal(n,sa,sx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to 1.
c     jack dongarra, linpack, 3/11/78.
c     modified to correct problem with negative increment, 8/21/90.
c
      real sa,sx(1)
      integer i,incx,ix,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      do 10 i = 1,n
        sx(ix) = sa*sx(ix)
        ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sx(i) = sa*sx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        sx(i) = sa*sx(i)
        sx(i + 1) = sa*sx(i + 1)
        sx(i + 2) = sa*sx(i + 2)
        sx(i + 3) = sa*sx(i + 3)
        sx(i + 4) = sa*sx(i + 4)
   50 continue
      return
      end
      integer function isamax(n,sx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c     modified to correct problem with negative increment, 8/21/90.
c
      real sx(1),smax
      integer i,incx,ix,n
c
      isamax = 0
      if( n .lt. 1 ) return
      isamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      smax = abs(sx(ix))
      ix = ix + incx
      do 10 i = 2,n
         if(abs(sx(ix)).le.smax) go to 5
         isamax = i
         smax = abs(sx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 smax = abs(sx(1))
      do 30 i = 2,n
         if(abs(sx(i)).le.smax) go to 30
         isamax = i
         smax = abs(sx(i))
   30 continue
      return
      end
      subroutine dgesl(a,lda,n,ipvt,b,job)
      integer lda,n,ipvt(1),job
      double precision a(lda,1),b(1)
c
c     dgesl solves the double precision system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgeco or dgefa.
c
c     on entry
c
c        a       double precision(lda, n)
c                the output from dgeco or dgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from dgeco or dgefa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgeco has set rcond .gt. 0.0
c        or dgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c
c     internal variables
c
      double precision ddot,t
      integer k,kb,l,nm1
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call daxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            t = ddot(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end
      double precision function ddot(n,dx,incx,dy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
      subroutine sgesl(a,lda,n,ipvt,b,job)
      integer lda,n,ipvt(1),job
      real a(lda,1),b(1)
c
c     sgesl solves the real system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by sgeco or sgefa.
c
c     on entry
c
c        a       real(lda, n)
c                the output from sgeco or sgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from sgeco or sgefa.
c
c        b       real(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if sgeco has set rcond .gt. 0.0
c        or sgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call sgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call sgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas saxpy,sdot
c
c     internal variables
c
      real sdot,t
      integer k,kb,l,nm1
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call saxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call saxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            t = sdot(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end
      real function sdot(n,sx,incx,sy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),stemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      stemp = 0.0e0
      sdot = 0.0e0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        stemp = stemp + sx(ix)*sy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      sdot = stemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        stemp = stemp + sx(i)*sy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +
     *   sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
   50 continue
   60 sdot = stemp
      return
      end
C-END-OF-FILE
cat > poiss.f <<C-END-OF-FILE
C  FACR(0) HELMHOLTZ-SOLVER FOR RECTANGLE WITH
C  DIRICHLET BOUNDARY CONDITIONS
C
      SUBROUTINE POISS(X,RADIUS,RHO,WORK,TRIGS,IFAX,NX,MX)
      DIMENSION X(*),WORK(*),TRIGS(*),IFAX(*)
C
      N=NX-1
      M=MX-1
      NL=N-1
      ML=M-1
C
C  SINE TRANSFORMS
C  ---------------
      CALL FFT44(X(NX+1),WORK,TRIGS,IFAX,1,NX,N,ML,-1)
C
C  COMPUTE COEFFICIENTS
C  --------------------
      DEL=2.0*ASIN(1.0)/FLOAT(N)
      RHO2=RHO*RHO
      DO 10 I=1,NL
      WORK(I+1)=2.*RHO2*COS(FLOAT(I)*DEL)-2.*(1.+RHO2)-RADIUS
   10 CONTINUE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 15 I=1,NL
      WORK(NX+I+1)=1.0/WORK(I+1)
   15 CONTINUE
      IA=NX+1
      DO 30 J=2,ML
      IB=IA+NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 20 I=1,NL
      WORK(IB+I)=1.0/(WORK(I+1)-WORK(IA+I))
   20 CONTINUE
      IA=IB
   30 CONTINUE
C
C  ZERO BOUNDARIES
C  ---------------
      LA=M*NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 40 I=1,NX
      X(I)=0.0
      X(LA+I)=0.0
   40 CONTINUE
      DO 50 I=NX,LA,NX
      X(I)=0.0
   50 CONTINUE
C
C  SOLVE TRIDIAGONAL SYSTEMS
C  -------------------------
      IA=1
      DO 70 J=1,ML
      IB=IA+NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 60 I=1,NL
      X(IB+I)=WORK(IB+I)*(X(IB+I)-X(IA+I))
   60 CONTINUE
      IA=IB
   70 CONTINUE
      IA=LA+1
      DO 90 J=1,ML
      IB=IA-NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 80 I=1,NL
      X(IB+I)=X(IB+I)-WORK(IB+I)*X(IA+I)
   80 CONTINUE
      IA=IB
   90 CONTINUE
C
C  INVERSE SINE TRANSFORMS
C  -----------------------
      CALL FFT44(X(NX+1),WORK,TRIGS,IFAX,1,NX,N,ML,+1)
C
      RETURN
      END
C
C  FACR(0) HELMHOLTZ-SOLVER IN SPHERICAL COORDINATES WITH
C  DIRICHLET BOUNDARY CONDITIONS.
C
      SUBROUTINE POISSA(X,RADIUS,WORK,TRIGS,IFAX,NX,MX,
     #                  C1,C2,C3)
      DIMENSION X(NX,MX),WORK(NX,MX),TRIGS(*),IFAX(*)
      DIMENSION C1(MX),C2(MX),C3(MX)
C
      N=NX-1
      M=MX-1
      ML=M-1
      DEL=2.0*ASIN(1.0)/FLOAT(N)
      PI=4.*ATAN(1.)
C
C  SINE TRANSFORMS
C  ---------------
      CALL FFT44(X(1,2),WORK,TRIGS,IFAX,1,NX,N,ML,-1)
C
C  COMPUTE COEFFICIENTS
C  --------------------

      DO 5 I=2,N
         WORK(I,1)=0.
         WORK(I,MX)=0.
5     CONTINUE

      DO 25 J=2,M
CDIR$ IVDEP
CVD$ NOSYNC
      DO 10 I=2,N
       AK=FLOAT(I-1)/FLOAT(N)
       B=(C2(J)+C1(J)+C3(J)*(2.-2.*COS(AK*PI))+RADIUS)
       WORK(I,J)=C1(J)/(B-C2(J)*WORK(I,J-1))
10    CONTINUE
25    CONTINUE

C  ZERO BOUNDARIES
C  ---------------

      DO 40 I=1,NX
      X(I,1)=0.0
      X(I,MX)=0.0
   40 CONTINUE
      DO 50 J=2,M
      X(1,J)=0.0
      X(NX,J)=0.0
   50 CONTINUE

C  SOLVE TRIDIAGONAL SYSTEMS
C  -------------------------

      DO 70 J=2,M
CDIR$ IVDEP
CVD$ NOSYNC
      DO 60 I=2,N
         X(I,J)=(WORK(I,J)/C1(J))*(-X(I,J)+C2(J)*X(I,J-1))
60    CONTINUE
70    CONTINUE

       DO 90 J=1,ML
         JJ=MX-J
CDIR$ IVDEP
CVD$ NOSYNC
        DO 80 I=2,N
          X(I,JJ)=X(I,JJ)+WORK(I,JJ)*X(I,JJ+1)
80      CONTINUE
90     CONTINUE

C  INVERSE SINE TRANSFORMS
C  -----------------------
500      CALL FFT44(X(1,2),WORK,TRIGS,IFAX,1,NX,N,ML,+1)

      RETURN
      END
C
      SUBROUTINE FAX(IFAX,N,MODE)
      DIMENSION IFAX(10)
      NN=N
      IF (IABS(MODE).EQ.1) GO TO 10
      IF (IABS(MODE).EQ.8) GO TO 10
      NN=N/2
      IF ((NN+NN).EQ.N) GO TO 10
      IFAX(1)=-99
      RETURN
   10 K=1
C     TEST FOR FACTORS OF 4
   20 IF (MOD(NN,4).NE.0) GO TO 30
      K=K+1
      IFAX(K)=4
      NN=NN/4
      IF (NN.EQ.1) GO TO 80
      GO TO 20
C     TEST FOR EXTRA FACTOR OF 2
   30 IF (MOD(NN,2).NE.0) GO TO 40
      K=K+1
      IFAX(K)=2
      NN=NN/2
      IF (NN.EQ.1) GO TO 80
C     TEST FOR FACTORS OF 3
   40 IF (MOD(NN,3).NE.0) GO TO 50
      K=K+1
      IFAX(K)=3
      NN=NN/3
      IF (NN.EQ.1) GO TO 80
      GO TO 40
C     NOW FIND REMAINING FACTORS
   50 L=5
      INC=2
C     INC ALTERNATELY TAKES ON VALUES 2 AND 4
   60 IF (MOD(NN,L).NE.0) GO TO 70
      K=K+1
      IFAX(K)=L
      NN=NN/L
      IF (NN.EQ.1) GO TO 80
      GO TO 60
   70 L=L+INC
      INC=6-INC
      GO TO 60
   80 IFAX(1)=K-1
C     IFAX(1) CONTAINS NUMBER OF FACTORS
C     IFAX(1) CONTAINS NUMBER OF FACTORS
      NFAX=IFAX(1)
C     SORT FACTORS INTO ASCENDING ORDER
      IF (NFAX.EQ.1) GO TO 110
      DO 100 II=2,NFAX
      ISTOP=NFAX+2-II
      DO 90 I=2,ISTOP
      IF (IFAX(I+1).GE.IFAX(I)) GO TO 90
      ITEM=IFAX(I)
      IFAX(I)=IFAX(I+1)
      IFAX(I+1)=ITEM
   90 CONTINUE
  100 CONTINUE
  110 CONTINUE
      RETURN
      END
C
      SUBROUTINE FFTRIG(TRIGS,N,MODE)
      DIMENSION TRIGS(*)
      PI=2.0*ASIN(1.0)
      IMODE=IABS(MODE)
      NN=N
      IF (IMODE.GT.1.AND.IMODE.LT.6) NN=N/2
      DEL=(PI+PI)/FLOAT(NN)
      L=NN+NN
      DO 10 I=1,L,2
      ANGLE=0.5*FLOAT(I-1)*DEL
      TRIGS(I)=COS(ANGLE)
      TRIGS(I+1)=SIN(ANGLE)
   10 CONTINUE
      IF (IMODE.EQ.1) RETURN
      IF (IMODE.EQ.8) RETURN
      DEL=0.5*DEL
      NH=(NN+1)/2
      L=NH+NH
      LA=NN+NN
      DO 20 I=1,L,2
      ANGLE=0.5*FLOAT(I-1)*DEL
      TRIGS(LA+I)=COS(ANGLE)
      TRIGS(LA+I+1)=SIN(ANGLE)
   20 CONTINUE
      IF (IMODE.LE.3) RETURN
      DEL=0.5*DEL
      LA=LA+NN
      IF (MODE.EQ.5) GO TO 40
      DO 30 I=2,NN
      ANGLE=FLOAT(I-1)*DEL
      TRIGS(LA+I)=2.0*SIN(ANGLE)
   30 CONTINUE
      RETURN
   40 CONTINUE
      DEL=0.5*DEL
      DO 50 I=2,N
      ANGLE=FLOAT(I-1)*DEL
      TRIGS(LA+I)=SIN(ANGLE)
   50 CONTINUE
      RETURN
      END
C
C     SUBROUTINE "FFT44" - MULTIPLE FAST SINE TRANSFORM
C     PROCEDURE USED TO CONVERT TO CONVERT TO HALF-LENGTH COMPLEX
C     TRANSFORMS IS THE INVERSE OF THAT GIVEN BY COOLEY, LEWIS ' WELCH
C     (J. SOUND VIB., VOL. 12 (1970), 315-337)
C
C     A IS THE ARRAY CONTAINING INPUT ' OUTPUT DATA
C     WORK IS AN AREA OF SIZE (N+1)*LOT
C     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
C     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2
C     INC IS THE INCREMENT WITHIN EACH DATA "VECTOR"
C         (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
C     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
C     N IS THE LENGTH OF THE DATA VECTORS (INCLUDING A(1)=0.0)
C     LOT IS THE NUMBER OF DATA VECTORS
C     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
C           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
C
C     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
C         PARALLEL
C
C     *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER
C
C     DEFINITION OF TRANSFORMS:
C     -------------------------
C
C     ISIGN=+1: X(J)=SUM(K=0,1,...,N-1)(B(K)*SIN(J*K*PI/N)),
C                                            J=0,1,...,N-1
C
C     ISIGN=-1: B(K)=(2/N)*SUM(J=0,1,...,N-1)(X(J)*SIN(J*K*PI/N)))
C                                            K=0,1,...,N-1
C
      SUBROUTINE FFT44(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
      DIMENSION A(*),WORK(*),TRIGS(*),IFAX(*)
C
      NH=N/2
      NH3=3*N/2
      NX=N+1
      INK=INC+INC
C
C     PREPROCESSING
C     -------------
C
      SCALE=0.25
      IF (ISIGN.EQ.-1) SCALE=0.5/FLOAT(N)
C
      IABASE=INC+1
      IBBASE=(N-1)*INC+1
      JABASE=2
      JBBASE=N
C
      DO 20 I=2,NH
      IA=IABASE
      IB=IBBASE
      JA=JABASE
      JB=JBBASE
      SINK=TRIGS(NH3+I)
CDIR$ IVDEP
CVD$ NOSYNC
      DO 10 J=1,LOT
      WORK(JA)=SCALE*(SINK*(A(IA)+A(IB))+(A(IA)-A(IB)))
      WORK(JB)=SCALE*(SINK*(A(IA)+A(IB))-(A(IA)-A(IB)))
      IA=IA+JUMP
      IB=IB+JUMP
      JA=JA+NX
      JB=JB+NX
   10 CONTINUE
      IABASE=IABASE+INC
      IBBASE=IBBASE-INC
      JABASE=JABASE+1
      JBBASE=JBBASE-1
   20 CONTINUE
      IA=IABASE
      JA=1
      JB=NH+1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 30 J=1,LOT
      WORK(JA)=0.0
      WORK(JB)=4.0*SCALE*A(IA)
      IA=IA+JUMP
      JA=JA+NX
      JB=JB+NX
   30 CONTINUE
C
C     COMPLEX FFT
C     -----------
C
      NFAX=IFAX(1)
      LA=1
      IGO=40
      DO 70 K=1,NFAX
      IF (IGO.EQ.50) GO TO 50
   40 CONTINUE
      CALL VPASSM(WORK(1),WORK(2),A(1),A(INC+1),TRIGS,2,INK,NX,JUMP,LOT,
     *    NH,IFAX(K+1),LA)
      IGO=50
      GO TO 60
   50 CONTINUE
      CALL VPASSM(A(1),A(INC+1),WORK(1),WORK(2),TRIGS,INK,2,JUMP,NX,LOT,
     *    NH,IFAX(K+1),LA)
      IGO=40
   60 CONTINUE
      LA=LA*IFAX(K+1)
   70 CONTINUE
C
C     POSTPROCESSING - FIRST STAGE
C     ----------------------------
C
      IF (MOD(NFAX,2).EQ.1) GO TO 80
      CALL FFT44A(WORK(1),WORK(2),A(1),A(INC+1),TRIGS,2,INK,NX,JUMP,
     *    LOT,N)
      GO TO 90
   80 CONTINUE
      CALL FFT44A(A(1),A(INC+1),WORK(1),WORK(2),TRIGS,INK,2,JUMP,NX,
     *    LOT,N)
C
C    POSTPROCESSING - SECOND STAGE
C    -----------------------------
C
   90 CONTINUE
      IF (MOD(NFAX,2).EQ.1) GO TO 100
      CALL FFT44B(A,A,INC,INC,JUMP,JUMP,LOT,N,0)
      GO TO 110
  100 CONTINUE
      CALL FFT44B(WORK,A,1,INC,NX,JUMP,LOT,N,1)
  110 CONTINUE
C
  120 FORMAT(1H0,8E14.5)
      RETURN
      END
C     SUBROUTINE FFT44A - FIRST STAGE OF POSTPROCESSING FOR
C     MULTIPLE SINE TRANSFORM
C
      SUBROUTINE FFT44A(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N)
      DIMENSION A(*),B(*),C(*),D(*),TRIGS(*)
C
      IA=1
      JA=1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 10 J=1,LOT
      C(JA)=A(IA)-B(IA)
      D(JA)=A(IA)+B(IA)
      IA=IA+INC3
      JA=JA+INC4
   10 CONTINUE
C
      NH=N/2
      IF (MOD(NH,2).EQ.1) GO TO 30
C
      IB=(NH*INC1)/2+1
      JB=(NH*INC2)/2+1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 20 J=1,LOT
      C(JB)=2.0*B(IB)
      D(JB)=2.0*A(IB)
      IB=IB+INC3
      JB=JB+INC4
   20 CONTINUE
C
   30 CONTINUE
      IABASE=INC1+1
      IBBASE=(NH-1)*INC1+1
      JABASE=INC2+1
      JBBASE=(NH-1)*INC2+1
      DO 50 K=3,NH,2
      IA=IABASE
      IB=IBBASE
      JA=JABASE
      JB=JBBASE
      CO=TRIGS(N+K)
      SI=TRIGS(N+K+1)
CDIR$ IVDEP
CVD$ NOSYNC
      DO 40 J=1,LOT
      C(JA)=(SI*(B(IA)+B(IB))-CO*(A(IA)-A(IB)))+(B(IA)-B(IB))
      C(JB)=(SI*(B(IA)+B(IB))-CO*(A(IA)-A(IB)))-(B(IA)-B(IB))
      D(JA)=(A(IA)+A(IB))+(SI*(A(IA)-A(IB))+CO*(B(IA)+B(IB)))
      D(JB)=(A(IA)+A(IB))-(SI*(A(IA)-A(IB))+CO*(B(IA)+B(IB)))
      IA=IA+INC3
      IB=IB+INC3
      JA=JA+INC4
      JB=JB+INC4
   40 CONTINUE
      IABASE=IABASE+INC1
      IBBASE=IBBASE-INC1
      JABASE=JABASE+INC2
      JBBASE=JBBASE-INC2
   50 CONTINUE
C
      RETURN
      END
C     SUBROUTINE FFT44B - SECOND STAGE OF POSTPROCESSING FOR
C     MULTIPLE SINE TRANSFORM   (ISEP=0 IF A,B ARE THE SAME ARRAY)
C
      SUBROUTINE FFT44B(A,B,INC1,INC2,INC3,INC4,LOT,N,ISEP)
      DIMENSION A(*),B(*)
C
      IA=1
      IB=(N-2)*INC1+1
      JA=1
      JB=(N-2)*INC2+1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 10 J=1,LOT
      B(JB+INC2)=-A(IA)
      B(JA+INC2)=A(IA+INC1)
      B(JA)=0.0
      IA=IA+INC3
      IB=IB+INC3
      JA=JA+INC4
      JB=JB+INC4
   10 CONTINUE
C
      IBASE=3*INC1+1
      JBASE=3*INC2+1
      INK1=INC1+INC1
      INK2=INC2+INC2
      NH=N/2
      DO 30 K=3,NH
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 20 L=1,LOT
      B(J)=B(J-INK2)+A(I)
      I=I+INC3
      J=J+INC4
   20 CONTINUE
      IBASE=IBASE+INK1
      JBASE=JBASE+INK2
   30 CONTINUE
C
      IF (ISEP.EQ.0) GO TO 60
C
      IBASE=2*INC1+1
      JBASE=2*INC2+1
      DO 50 K=2,NH
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 40 L=1,LOT
      B(J)=A(I)
      I=I+INC3
      J=J+INC4
   40 CONTINUE
      IBASE=IBASE+INK1
      JBASE=JBASE+INK2
   50 CONTINUE
C
   60 CONTINUE
      RETURN
      END
C
C     SUBROUTINE 'VPASSM' - MULTIPLE VERSION OF 'VPASSA'
C     PERFORMS ONE PASS THROUGH DATA
C     AS PART OF MULTIPLE COMPLEX FFT ROUTINE
C     A IS FIRST REAL INPUT VECTOR
C     B IS FIRST IMAGINARY INPUT VECTOR
C     C IS FIRST REAL OUTPUT VECTOR
C     D IS FIRST IMAGINARY OUTPUT VECTOR
C     TRIGS IS PRECALCULATED TABLE OF SINES ' COSINES
C     INC1 IS ADDRESSING INCREMENT FOR A AND B
C     INC2 IS ADDRESSING INCREMENT FOR C AND D
C     INC3 IS ADDRESSING INCREMENT BETWEEN A'S & B'S
C     INC4 IS ADDRESSING INCREMENT BETWEEN C'S & D'S
C     LOT IS THE NUMBER OF VECTORS
C     N IS LENGTH OF VECTORS
C     IFAC IS CURRENT FACTOR OF N
C     LA IS PRODUCT OF PREVIOUS FACTORS
C
      SUBROUTINE VPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA)
      DIMENSION A(*),B(*),C(*),D(*),TRIGS(*)
      DATA SIN36/0.587785252292473/,COS36/0.809016994374947/,
     *     SIN72/0.951056516295154/,COS72/0.309016994374947/,
     *     SIN60/0.866025403784437/
C
      M=N/IFAC
      IINK=M*INC1
      JINK=LA*INC2
      JUMP=(IFAC-1)*JINK
      IBASE=0
      JBASE=0
      IGO=IFAC-1
      IF (IGO.GT.4) RETURN
      GO TO (10,50,90,130),IGO
C
C     CODING FOR FACTOR 2
C
   10 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      DO 20 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 15 IJK=1,LOT
      C(JA+J)=A(IA+I)+A(IB+I)
      D(JA+J)=B(IA+I)+B(IB+I)
      C(JB+J)=A(IA+I)-A(IB+I)
      D(JB+J)=B(IA+I)-B(IB+I)
      I=I+INC3
      J=J+INC4
   15 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   20 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 40 K=LA1,M,LA
      KB=K+K-2
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      DO 30 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 25 IJK=1,LOT
      C(JA+J)=A(IA+I)+A(IB+I)
      D(JA+J)=B(IA+I)+B(IB+I)
      C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I))
      D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I))
      I=I+INC3
      J=J+INC4
   25 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   30 CONTINUE
      JBASE=JBASE+JUMP
   40 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 3
C
   50 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      DO 60 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 55 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
      C(JB+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))
      C(JC+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))
      D(JB+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))
      D(JC+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))
      I=I+INC3
      J=J+INC4
   55 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   60 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 80 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      DO 70 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 65 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
      C(JB+J)=
     *    C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   -S1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      D(JB+J)=
     *    S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   +C1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      C(JC+J)=
     *    C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   -S2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      D(JC+J)=
     *    S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   +C2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      I=I+INC3
      J=J+INC4
   65 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   70 CONTINUE
      JBASE=JBASE+JUMP
   80 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 4
C
   90 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      ID=IC+IINK
      JD=JC+JINK
      DO 100 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 95 IJK=1,LOT
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
      C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
      D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))
      C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))
      C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))
      D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))
      D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))
      I=I+INC3
      J=J+INC4
   95 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  100 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 120 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      KD=KC+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      C3=TRIGS(KD+1)
      S3=TRIGS(KD+2)
      DO 110 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 105 IJK=1,LOT
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
      C(JC+J)=
     *    C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
     *   -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
      D(JC+J)=
     *    S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
     *   +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
      C(JB+J)=
     *    C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))
     *   -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
      D(JB+J)=
     *    S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))
     *   +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
      C(JD+J)=
     *    C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))
     *   -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
      D(JD+J)=
     *    S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))
     *   +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
      I=I+INC3
      J=J+INC4
  105 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  110 CONTINUE
      JBASE=JBASE+JUMP
  120 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 5
C
  130 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      ID=IC+IINK
      JD=JC+JINK
      IE=ID+IINK
      JE=JD+JINK
      DO 140 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 135 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
      C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *  -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
      C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *  +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
      D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *  +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
      D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *  -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
      C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *  -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
      C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *  +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
      D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *  +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
      D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *  -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
      I=I+INC3
      J=J+INC4
  135 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  140 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 160 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      KD=KC+KB
      KE=KD+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      C3=TRIGS(KD+1)
      S3=TRIGS(KD+2)
      C4=TRIGS(KE+1)
      S4=TRIGS(KE+2)
      DO 150 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 145 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
      C(JB+J)=
     *    C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      D(JB+J)=
     *    S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      C(JE+J)=
     *    C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      D(JE+J)=
     *    S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      C(JC+J)=
     *    C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      D(JC+J)=
     *    S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      C(JD+J)=
     *    C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      D(JD+J)=
     *    S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      I=I+INC3
      J=J+INC4
  145 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  150 CONTINUE
      JBASE=JBASE+JUMP
  160 CONTINUE
      RETURN
      END
C
C  FACR(0) POISSON-SOLVER FOR RECTANGLE WITH
C  DIRICHLET BOUNDARY CONDITIONS
C
      SUBROUTINE DPOISS(X,RADIUS,RHO,WORK,TRIGS,IFAX,NX,MX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(*),WORK(*),TRIGS(*),IFAX(*)
C
      N=NX-1
      M=MX-1
      NL=N-1
      ML=M-1
C
C  SINE TRANSFORMS
C  ---------------
      CALL DFFT44(X(NX+1),WORK,TRIGS,IFAX,1,NX,N,ML,-1)
C
C  COMPUTE COEFFICIENTS
C  --------------------
      DEL=2.0D0*ASIN(1.0D0)/FLOAT(N)
      RHO2=RHO*RHO
      DO 10 I=1,NL
      WORK(I+1)=2.D0*RHO2*COS(FLOAT(I)*DEL)-2.D0*(1.D0+RHO2)-RADIUS
   10 CONTINUE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 15 I=1,NL
      WORK(NX+I+1)=1.0D0/WORK(I+1)
   15 CONTINUE
      IA=NX+1
      DO 30 J=2,ML
      IB=IA+NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 20 I=1,NL
      WORK(IB+I)=1.0D0/(WORK(I+1)-WORK(IA+I))
   20 CONTINUE
      IA=IB
   30 CONTINUE
C
C  ZERO BOUNDARIES
C  ---------------
      LA=M*NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 40 I=1,NX
      X(I)=0.0D0
      X(LA+I)=0.0D0
   40 CONTINUE
      DO 50 I=NX,LA,NX
      X(I)=0.0D0
   50 CONTINUE
C
C  SOLVE TRIDIAGONAL SYSTEMS
C  -------------------------
      IA=1
      DO 70 J=1,ML
      IB=IA+NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 60 I=1,NL
      X(IB+I)=WORK(IB+I)*(X(IB+I)-X(IA+I))
   60 CONTINUE
      IA=IB
   70 CONTINUE
      IA=LA+1
      DO 90 J=1,ML
      IB=IA-NX
CDIR$ IVDEP
CVD$ NOSYNC
      DO 80 I=1,NL
      X(IB+I)=X(IB+I)-WORK(IB+I)*X(IA+I)
   80 CONTINUE
      IA=IB
   90 CONTINUE
C
C  INVERSE SINE TRANSFORMS
C  -----------------------
      CALL DFFT44(X(NX+1),WORK,TRIGS,IFAX,1,NX,N,ML,+1)
C
      RETURN
      END
C
      SUBROUTINE DFAX(IFAX,N,MODE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION IFAX(10)
      NN=N
      IF (IABS(MODE).EQ.1) GO TO 10
      IF (IABS(MODE).EQ.8) GO TO 10
      NN=N/2
      IF ((NN+NN).EQ.N) GO TO 10
      IFAX(1)=-99
      RETURN
   10 K=1
C     TEST FOR FACTORS OF 4
   20 IF (MOD(NN,4).NE.0) GO TO 30
      K=K+1
      IFAX(K)=4
      NN=NN/4
      IF (NN.EQ.1) GO TO 80
      GO TO 20
C     TEST FOR EXTRA FACTOR OF 2
   30 IF (MOD(NN,2).NE.0) GO TO 40
      K=K+1
      IFAX(K)=2
      NN=NN/2
      IF (NN.EQ.1) GO TO 80
C     TEST FOR FACTORS OF 3
   40 IF (MOD(NN,3).NE.0) GO TO 50
      K=K+1
      IFAX(K)=3
      NN=NN/3
      IF (NN.EQ.1) GO TO 80
      GO TO 40
C     NOW FIND REMAINING FACTORS
   50 L=5
      INC=2
C     INC ALTERNATELY TAKES ON VALUES 2 AND 4
   60 IF (MOD(NN,L).NE.0) GO TO 70
      K=K+1
      IFAX(K)=L
      NN=NN/L
      IF (NN.EQ.1) GO TO 80
      GO TO 60
   70 L=L+INC
      INC=6-INC
      GO TO 60
   80 IFAX(1)=K-1
C     IFAX(1) CONTAINS NUMBER OF FACTORS
C     IFAX(1) CONTAINS NUMBER OF FACTORS
      NFAX=IFAX(1)
C     SORT FACTORS INTO ASCENDING ORDER
      IF (NFAX.EQ.1) GO TO 110
      DO 100 II=2,NFAX
      ISTOP=NFAX+2-II
      DO 90 I=2,ISTOP
      IF (IFAX(I+1).GE.IFAX(I)) GO TO 90
      ITEM=IFAX(I)
      IFAX(I)=IFAX(I+1)
      IFAX(I+1)=ITEM
   90 CONTINUE
  100 CONTINUE
  110 CONTINUE
      RETURN
      END
C
      SUBROUTINE DFFTRG(TRIGS,N,MODE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TRIGS(*)
      PI=2.0D0*ASIN(1.0D0)
      IMODE=IABS(MODE)
      NN=N
      IF (IMODE.GT.1.AND.IMODE.LT.6) NN=N/2
      DEL=(PI+PI)/FLOAT(NN)
      L=NN+NN
      DO 10 I=1,L,2
      ANGLE=0.5D0*FLOAT(I-1)*DEL
      TRIGS(I)=COS(ANGLE)
      TRIGS(I+1)=SIN(ANGLE)
   10 CONTINUE
      IF (IMODE.EQ.1) RETURN
      IF (IMODE.EQ.8) RETURN
      DEL=0.5D0*DEL
      NH=(NN+1)/2
      L=NH+NH
      LA=NN+NN
      DO 20 I=1,L,2
      ANGLE=0.5D0*FLOAT(I-1)*DEL
      TRIGS(LA+I)=COS(ANGLE)
      TRIGS(LA+I+1)=SIN(ANGLE)
   20 CONTINUE
      IF (IMODE.LE.3) RETURN
      DEL=0.5D0*DEL
      LA=LA+NN
      IF (MODE.EQ.5) GO TO 40
      DO 30 I=2,NN
      ANGLE=FLOAT(I-1)*DEL
      TRIGS(LA+I)=2.0D0*SIN(ANGLE)
   30 CONTINUE
      RETURN
   40 CONTINUE
      DEL=0.5D0*DEL
      DO 50 I=2,N
      ANGLE=FLOAT(I-1)*DEL
      TRIGS(LA+I)=SIN(ANGLE)
   50 CONTINUE
      RETURN
      END
C
C     SUBROUTINE "DFFT44" - MULTIPLE FAST SINE TRANSFORM
C     PROCEDURE USED TO CONVERT TO CONVERT TO HALF-LENGTH COMPLEX
C     TRANSFORMS IS THE INVERSE OF THAT GIVEN BY COOLEY, LEWIS ' WELCH
C     (J. SOUND VIB., VOL. 12 (1970), 315-337)
C
C     A IS THE ARRAY CONTAINING INPUT ' OUTPUT DATA
C     WORK IS AN AREA OF SIZE (N+1)*LOT
C     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
C     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2
C     INC IS THE INCREMENT WITHIN EACH DATA "VECTOR"
C         (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
C     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
C     N IS THE LENGTH OF THE DATA VECTORS (INCLUDING A(1)=0.0)
C     LOT IS THE NUMBER OF DATA VECTORS
C     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
C           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
C
C     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
C         PARALLEL
C
C     *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER
C
C     DEFINITION OF TRANSFORMS:
C     -------------------------
C
C     ISIGN=+1: X(J)=SUM(K=0,1,...,N-1)(B(K)*SIN(J*K*PI/N)),
C                                            J=0,1,...,N-1
C
C     ISIGN=-1: B(K)=(2/N)*SUM(J=0,1,...,N-1)(X(J)*SIN(J*K*PI/N)))
C                                            K=0,1,...,N-1
C
      SUBROUTINE DFFT44(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(*),WORK(*),TRIGS(*),IFAX(*)
C
      NH=N/2
      NH3=3*N/2
      NX=N+1
      INK=INC+INC
C
C     PREPROCESSING
C     -------------
C
      SCALE=0.25D0
      IF (ISIGN.EQ.-1) SCALE=0.5D0/FLOAT(N)
C
      IABASE=INC+1
      IBBASE=(N-1)*INC+1
      JABASE=2
      JBBASE=N
C
      DO 20 I=2,NH
      IA=IABASE
      IB=IBBASE
      JA=JABASE
      JB=JBBASE
      SINK=TRIGS(NH3+I)
CDIR$ IVDEP
CVD$ NOSYNC
      DO 10 J=1,LOT
      WORK(JA)=SCALE*(SINK*(A(IA)+A(IB))+(A(IA)-A(IB)))
      WORK(JB)=SCALE*(SINK*(A(IA)+A(IB))-(A(IA)-A(IB)))
      IA=IA+JUMP
      IB=IB+JUMP
      JA=JA+NX
      JB=JB+NX
   10 CONTINUE
      IABASE=IABASE+INC
      IBBASE=IBBASE-INC
      JABASE=JABASE+1
      JBBASE=JBBASE-1
   20 CONTINUE
      IA=IABASE
      JA=1
      JB=NH+1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 30 J=1,LOT
      WORK(JA)=0.0D0
      WORK(JB)=4.0D0*SCALE*A(IA)
      IA=IA+JUMP
      JA=JA+NX
      JB=JB+NX
   30 CONTINUE
C
C     COMPLEX FFT
C     -----------
C
      NFAX=IFAX(1)
      LA=1
      IGO=40
      DO 70 K=1,NFAX
      IF (IGO.EQ.50) GO TO 50
   40 CONTINUE
      CALL DVPASM(WORK(1),WORK(2),A(1),A(INC+1),TRIGS,2,INK,NX,JUMP,
     *    LOT,NH,IFAX(K+1),LA)
      IGO=50
      GO TO 60
   50 CONTINUE
      CALL DVPASM(A(1),A(INC+1),WORK(1),WORK(2),TRIGS,INK,2,JUMP,NX,
     *    LOT,NH,IFAX(K+1),LA)
      IGO=40
   60 CONTINUE
      LA=LA*IFAX(K+1)
   70 CONTINUE
C
C     POSTPROCESSING - FIRST STAGE
C     ----------------------------
C
      IF (MOD(NFAX,2).EQ.1) GO TO 80
      CALL DFFT4A(WORK(1),WORK(2),A(1),A(INC+1),TRIGS,2,INK,NX,JUMP,
     *    LOT,N)
      GO TO 90
   80 CONTINUE
      CALL DFFT4A(A(1),A(INC+1),WORK(1),WORK(2),TRIGS,INK,2,JUMP,NX,
     *    LOT,N)
C
C    POSTPROCESSING - SECOND STAGE
C    -----------------------------
C
   90 CONTINUE
      IF (MOD(NFAX,2).EQ.1) GO TO 100
      CALL DFFT4B(A,A,INC,INC,JUMP,JUMP,LOT,N,0)
      GO TO 110
  100 CONTINUE
      CALL DFFT4B(WORK,A,1,INC,NX,JUMP,LOT,N,1)
  110 CONTINUE
C
  120 FORMAT(1H0,8D14.5)
      RETURN
      END
C     SUBROUTINE DFFT4A - FIRST STAGE OF POSTPROCESSING FOR
C     MULTIPLE SINE TRANSFORM
C
      SUBROUTINE DFFT4A(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(*),B(*),C(*),D(*),TRIGS(*)
C
      IA=1
      JA=1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 10 J=1,LOT
      C(JA)=A(IA)-B(IA)
      D(JA)=A(IA)+B(IA)
      IA=IA+INC3
      JA=JA+INC4
   10 CONTINUE
C
      NH=N/2
      IF (MOD(NH,2).EQ.1) GO TO 30
C
      IB=(NH*INC1)/2+1
      JB=(NH*INC2)/2+1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 20 J=1,LOT
      C(JB)=2.0D0*B(IB)
      D(JB)=2.0D0*A(IB)
      IB=IB+INC3
      JB=JB+INC4
   20 CONTINUE
C
   30 CONTINUE
      IABASE=INC1+1
      IBBASE=(NH-1)*INC1+1
      JABASE=INC2+1
      JBBASE=(NH-1)*INC2+1
      DO 50 K=3,NH,2
      IA=IABASE
      IB=IBBASE
      JA=JABASE
      JB=JBBASE
      CO=TRIGS(N+K)
      SI=TRIGS(N+K+1)
CDIR$ IVDEP
CVD$ NOSYNC
      DO 40 J=1,LOT
      C(JA)=(SI*(B(IA)+B(IB))-CO*(A(IA)-A(IB)))+(B(IA)-B(IB))
      C(JB)=(SI*(B(IA)+B(IB))-CO*(A(IA)-A(IB)))-(B(IA)-B(IB))
      D(JA)=(A(IA)+A(IB))+(SI*(A(IA)-A(IB))+CO*(B(IA)+B(IB)))
      D(JB)=(A(IA)+A(IB))-(SI*(A(IA)-A(IB))+CO*(B(IA)+B(IB)))
      IA=IA+INC3
      IB=IB+INC3
      JA=JA+INC4
      JB=JB+INC4
   40 CONTINUE
      IABASE=IABASE+INC1
      IBBASE=IBBASE-INC1
      JABASE=JABASE+INC2
      JBBASE=JBBASE-INC2
   50 CONTINUE
C
      RETURN
      END
C     SUBROUTINE DFFT4B - SECOND STAGE OF POSTPROCESSING FOR
C     MULTIPLE SINE TRANSFORM   (ISEP=0 IF A,B ARE THE SAME ARRAY)
C
      SUBROUTINE DFFT4B(A,B,INC1,INC2,INC3,INC4,LOT,N,ISEP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(*),B(*)
C
      IA=1
      IB=(N-2)*INC1+1
      JA=1
      JB=(N-2)*INC2+1
CDIR$ IVDEP
CVD$ NOSYNC
      DO 10 J=1,LOT
      B(JB+INC2)=-A(IA)
      B(JA+INC2)=A(IA+INC1)
      B(JA)=0.0D0
      IA=IA+INC3
      IB=IB+INC3
      JA=JA+INC4
      JB=JB+INC4
   10 CONTINUE
C
      IBASE=3*INC1+1
      JBASE=3*INC2+1
      INK1=INC1+INC1
      INK2=INC2+INC2
      NH=N/2
      DO 30 K=3,NH
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 20 L=1,LOT
      B(J)=B(J-INK2)+A(I)
      I=I+INC3
      J=J+INC4
   20 CONTINUE
      IBASE=IBASE+INK1
      JBASE=JBASE+INK2
   30 CONTINUE
C
      IF (ISEP.EQ.0) GO TO 60
C
      IBASE=2*INC1+1
      JBASE=2*INC2+1
      DO 50 K=2,NH
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 40 L=1,LOT
      B(J)=A(I)
      I=I+INC3
      J=J+INC4
   40 CONTINUE
      IBASE=IBASE+INK1
      JBASE=JBASE+INK2
   50 CONTINUE
C
   60 CONTINUE
      RETURN
      END
C
C     SUBROUTINE 'DVPASM' - MULTIPLE VERSION OF 'DVPASSA'
C     PERFORMS ONE PASS THROUGH DATA
C     AS PART OF MULTIPLE COMPLEX FFT ROUTINE
C     A IS FIRST REAL INPUT VECTOR
C     B IS FIRST IMAGINARY INPUT VECTOR
C     C IS FIRST REAL OUTPUT VECTOR
C     D IS FIRST IMAGINARY OUTPUT VECTOR
C     TRIGS IS PRECALCULATED TABLE OF SINES ' COSINES
C     INC1 IS ADDRESSING INCREMENT FOR A AND B
C     INC2 IS ADDRESSING INCREMENT FOR C AND D
C     INC3 IS ADDRESSING INCREMENT BETWEEN A'S & B'S
C     INC4 IS ADDRESSING INCREMENT BETWEEN C'S & D'S
C     LOT IS THE NUMBER OF VECTORS
C     N IS LENGTH OF VECTORS
C     IFAC IS CURRENT FACTOR OF N
C     LA IS PRODUCT OF PREVIOUS FACTORS
C
      SUBROUTINE DVPASM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,
     #                   IFAC,LA)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(*),B(*),C(*),D(*),TRIGS(*)
      DATA SIN36/0.587785252292473D0/,COS36/0.809016994374947D0/,
     *     SIN72/0.951056516295154D0/,COS72/0.309016994374947D0/,
     *     SIN60/0.866025403784437D0/
C
      M=N/IFAC
      IINK=M*INC1
      JINK=LA*INC2
      JUMP=(IFAC-1)*JINK
      IBASE=0
      JBASE=0
      IGO=IFAC-1
      IF (IGO.GT.4) RETURN
      GO TO (10,50,90,130),IGO
C
C     CODING FOR FACTOR 2
C
   10 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      DO 20 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 15 IJK=1,LOT
      C(JA+J)=A(IA+I)+A(IB+I)
      D(JA+J)=B(IA+I)+B(IB+I)
      C(JB+J)=A(IA+I)-A(IB+I)
      D(JB+J)=B(IA+I)-B(IB+I)
      I=I+INC3
      J=J+INC4
   15 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   20 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 40 K=LA1,M,LA
      KB=K+K-2
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      DO 30 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 25 IJK=1,LOT
      C(JA+J)=A(IA+I)+A(IB+I)
      D(JA+J)=B(IA+I)+B(IB+I)
      C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I))
      D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I))
      I=I+INC3
      J=J+INC4
   25 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   30 CONTINUE
      JBASE=JBASE+JUMP
   40 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 3
C
   50 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      DO 60 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 55 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
      C(JB+J)=
     # (A(IA+I)-0.5D0*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))
      C(JC+J)=
     # (A(IA+I)-0.5D0*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))
      D(JB+J)=
     # (B(IA+I)-0.5D0*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))
      D(JC+J)=
     # (B(IA+I)-0.5D0*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))
      I=I+INC3
      J=J+INC4
   55 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   60 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 80 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      DO 70 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 65 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
      C(JB+J)=
     * C1*((A(IA+I)-0.5D0*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *-S1*((B(IA+I)-0.5D0*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      D(JB+J)=
     *  S1*((A(IA+I)-0.5D0*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     * +C1*((B(IA+I)-0.5D0*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      C(JC+J)=
     *    C2*((A(IA+I)-0.5D0*(A(IB+I)+A(IC+I)))+
     #       (SIN60*(B(IB+I)-B(IC+I))))
     *   -S2*((B(IA+I)-0.5D0*(B(IB+I)+B(IC+I)))-
     #       (SIN60*(A(IB+I)-A(IC+I))))
      D(JC+J)=
     *    S2*((A(IA+I)-0.5D0*(A(IB+I)+A(IC+I)))+
     #       (SIN60*(B(IB+I)-B(IC+I))))
     *   +C2*((B(IA+I)-0.5D0*(B(IB+I)+B(IC+I)))-
     #       (SIN60*(A(IB+I)-A(IC+I))))
      I=I+INC3
      J=J+INC4
   65 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   70 CONTINUE
      JBASE=JBASE+JUMP
   80 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 4
C
   90 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      ID=IC+IINK
      JD=JC+JINK
      DO 100 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 95 IJK=1,LOT
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
      C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
      D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))
      C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))
      C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))
      D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))
      D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))
      I=I+INC3
      J=J+INC4
   95 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  100 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 120 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      KD=KC+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      C3=TRIGS(KD+1)
      S3=TRIGS(KD+2)
      DO 110 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 105 IJK=1,LOT
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
      C(JC+J)=
     *    C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
     *   -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
      D(JC+J)=
     *    S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
     *   +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
      C(JB+J)=
     *    C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))
     *   -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
      D(JB+J)=
     *    S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))
     *   +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
      C(JD+J)=
     *    C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))
     *   -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
      D(JD+J)=
     *    S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))
     *   +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
      I=I+INC3
      J=J+INC4
  105 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  110 CONTINUE
      JBASE=JBASE+JUMP
  120 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 5
C
  130 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      ID=IC+IINK
      JD=JC+JINK
      IE=ID+IINK
      JE=JD+JINK
      DO 140 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 135 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
      C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *  -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
      C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *  +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
      D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *  +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
      D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *  -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
      C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *  -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
      C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *  +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
      D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *  +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
      D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *  -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
      I=I+INC3
      J=J+INC4
  135 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  140 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 160 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      KD=KC+KB
      KE=KD+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      C3=TRIGS(KD+1)
      S3=TRIGS(KD+2)
      C4=TRIGS(KE+1)
      S4=TRIGS(KE+2)
      DO 150 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
CVD$ NOSYNC
      DO 145 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
      C(JB+J)=
     *    C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      D(JB+J)=
     *    S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      C(JE+J)=
     *    C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      D(JE+J)=
     *    S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      C(JC+J)=
     *    C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      D(JC+J)=
     *    S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      C(JD+J)=
     *    C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      D(JD+J)=
     *    S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      I=I+INC3
      J=J+INC4
  145 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  150 CONTINUE
      JBASE=JBASE+JUMP
  160 CONTINUE
      RETURN
      END
C-END-OF-FILE
cat > solvers.f <<C-END-OF-FILE
c***DECK : capc.f
      SUBROUTINE CAPC(F,U,WORK,TRIGS,NX,NY,HELM,DX,DY,
     #                C,CHARGE,IPS,VP,IP,JP,NP,ICFLAG,
     #                BA,BB,BC,BD,IERROR)
c
c     Programmer: P. Cummins
c     Version 1.0
c
c  *********************************************************************
c
c     CAPC is a routine to solve the second-order accurate finite
c     difference approximation to the Helmholtz equation
c
c     (d/dX)(dU/dX) + (d/dY)(dU/dY) - H*U = F(X,Y)
c
c     with Dirichlet boundary conditions in irregular polygonal
c     two-dimensional domains in Cartesian coordinates using the
c     capacitance matrix method.
c
c     ******************* ON ENTRY ********************
c
c         NX   integer constant
c           The number of mesh points in the X direction
c           Note that the value of NX is restricted such
c           that (NX-1) is of the form (2**I * 3**J * 5**K),
c           where I is an integer greater than or equal to
c           one, and J and K are integers greater than or
c           equal to zero. In addition, NX-1 must be greater
c           than or equal to 8.
c
c         NY   integer constant
c           The number of mesh points in the Y direction. NY-1
c           must be greater than or equal to 8.
c
c         NP   integer constant
c           The number of irregular boundary points. NP must be
c           greater than or equal to 1.
c
c          F   real array (NX,NY)
c            The source function on the right hand side. The value
c            of F is arbitrary at points in the rectangle which are
c            outside the irregular solution domain. This array is
c            unchanged on exit.
c
c       HELM   real constant
c            The coefficient H in the Helmholtz equation given above;
c            unchanged on exit.
c
c       WORK   real array (NX,NY)
c            A work space array
c
c     CHARGE   real array (NP)
c            A work space array
c
c         DX   real constant
c            The discretization interval in the X direction.
c            The condition DX > 0 must hold.
c
c         DY   real constant
c            The discretization interval in the Y direction.
c            The condition DY > 0 must hold.
c
c         BA   real array (NY)
c           The boundary values on the regular boundary
c           at (I=1, J=2,NY-1).
c
c         BB   real array (NY)
c           The boundary values on the regular boundary
c           at (I=NX, J=2,NY-1).
c
c         BC   real array (nx)
c           The boundary values on the regular boundary
c           at (I=1,NX, J=1).
c
c         BD   real array (NX)
c           The boundary values on the regular boundary
c           at (I=1,NX, J=NY).
c
c         VP   real array (NP)
c           The boundary values on the irregular boundary,
c           defined by the (IP,JP) coordinates.
c
c         IP   integer array (NP)
c           The coordinates in the X direction of the
c           irregular boundary points. The condition,
c           1 < IP(N) < NX, must hold for all N.
c
c         JP   integer array (NP)
c           The coordinates in the Y direction of the
c           irregular boundary points. The condition,
c           1 < JP(N) < NY, must hold for all N.
c
c            Additional notes on defining the IP and JP arrays are
c            given below.
c
c      Note: The order in which the irregular boundary grid points, defined
c            by the (IP(N),JP(N)) pairs, are specified is immaterial.
c            However, these grid points must all be distinct; there must
c            not be a repeated (IP(N),JP(N)) pair. Failure to ensure this
c            will result in the error condition IERROR=4, defined below.
c
c     ICFLAG   integer constant
c           A flag:  ICFLAG=0 - do preprocessing only and compute
c                               the matrix C and indices IPS.
c
c                    ICFLAG=1 - do preprocessing to compute the
c                               matrix C and indices IPS, and solve
c                               the Helmholtz problem.
c
c                    ICFLAG=2 - solve the Helmholtz problem only; this
c                               assumes that preprocessing has been
c                               done and that C and IPS are defined.
c
c            Any other value for ICFLAG gives IERROR=5 and no attempt
c            is made to do preprocessing or to solve. On return,
c            ICFLAG=2.
c
c
c     ******************* ON RETURN ********************
c
c          U   real array (NX,NY)
c           The solution in the irregular domain, with
c           U(IP(N),JP(N)) = VP(N) for N=1,NP on the irregular
c           portion of the boundary and also with
c              U(1,J)=BA(J) for J=2,NY-1
c              U(NX,J)=BB(J) for J=2,NY-1
c              U(I,1)=BC(I) for I=1,NX
c              U(I,NY)=BD(I) for I=1,NX.
c
c          C   real array (NP,NP)
c           The inverse capacitance matrix. This array must not
c           be altered on successive calls to the routine with
c           different source functions when ICFLAG=2.
c
c        IPS   integer array (NP)
c           An integer vector of pivot indices. This array must
c           not be alter on successive calls with different
c           source functions when ICFLAG=2.
c
c      TRIGS   real array (2*(NX-1))
c            An array of coefficients required by the FFT
c            routine employed by the fast solver. This array
c            is initialized on the first call to CAPC, or
c            whenever ICFLAG=0 or 1. The array must not be
c            altered between successive calls to the routine
c            with different source functions when ICFLAG=2.
c
c     IERROR   integer constant
c           An error flag on input variables. On return IERROR
c           has the following meaning:
c
c           IERROR=0  - normal status, no error detected.
c                 =1  - value of NX-1 or NY-1 is less than 8.
c                 =2  - illegal value for NX.
c                 =3  - illegal value for IP or JP.
c                 =4  - an error occurred in factoring the inverse
c                       capacitance matrix, i.e. the matrix is singular.
c                       This condition will occur if  the user has
c                       defined duplicate irregular boundary points.
c                 =5  - illegal value for ICFLAG.
c                 =6  - illegal value for DX or DY.
c
c     Note: Error checking is normally only done on the first call
c           to the routine. The user may force the routine to do error
c           checking at any time by setting IERROR=-1 on entry.
c
c
c   *********       NOTES ON CALLING ROUTINE CAPC       ********
c
c     In applying the capacitance matrix method to the solution
c     of the Helmholtz equation, routine CAPC works in two
c     stages. The first stage is a preprocessing one in which
c     the arrays C and IPS are determined. The second stage is
c     the actual solution of the problem for a given source
c     function. The flag ICFLAG determines the operation of the
c     routine:
c
c          If the routine is called with ICFLAG=0, only
c          preprocessing is done, the arrays C and IPC
c          are determined, and CAPC immediately returns to
c          the calling program without solving. On return
c          ICFLAG=2
c
c          If the routine is called with ICFLAG=1,
c          preprocessing is done, the arrays C and IPC
c          are determined, and the solution for the
c          given source function is determined. On return
c          ICFLAG=2.
c
c          If the routine is called with ICFLAG=2,
c          the solution for the given source function
c          is determined. This assumes that preprocessing
c          has been done and that the C and IPS arrays
c          are defined. On return ICFLAG=2.
c
c     If the user wishes to solve the Helmholtz equation in a
c     given domain repeatedly with different source functions,
c     then the routine could be called with ICFLAG=1 initially.
c     Afterwards this flag is set to 2 and need not be altered.
c     If the user has stored the arrays C and IPS from a previous
c     run, then these may be supplied as input data and the
c     routine may be called initially with ICFLAG=2 to obtain
c     the solution directly.
c
c     Preprocessing need only be done once for a given domain
c     geometry and Helmholtz coefficient. The user-specified
c     variables listed below must not change between successive
c     calls to CAPC with ICFLAG=2. If the user changes any
c     one of the following variables:
c
c                      NX,
c                      NY,
c                      DX,
c                      DY,
c                      IP,
c                      JP,
c                      NP,
c                      HELM,
c
c     then the preprocessing must be redone (ICFLAG=0 or 1). The
c     arrays TRIGS, C, and IPS should also remain unchanged on
c     successive calls with ICFLAG=2.
c
c
c   *********   NOTES ON DEFINING THE SOLUTION DOMAIN   *********
c
c   The routine solves an elliptic equation in an irregular,
c   polygonal domain. The domain is polygonal because its boundary
c   consists of straight line segments parallel or diagonal to
c   the discrete mesh. This domain is embedded in a rectangle
c   dimensioned NX by NY. The user is free to position the irregular
c   domain as desired within the embedding rectangle and the sides of
c   the rectangle may be part of the boundary of the solution domain.
c
c   The boundary of the irregular domain is described by a set of
c   grid points which are classified as 'regular' boundary points
c   and 'irregular' boundary points. A regular boundary point is one
c   which lies on one of the sides of the embedding rectangle. An
c   irregular boundary point is one which is not on a side of the
c   embedding rectangle. The routine requires that the user provide
c   the coordinates of the irregular boundary points through the
c   one-dimensional arrays IP and JP, both dimensioned NP. The user
c   must ensure that the specified coordinates of the irregular boundary
c   points do not lie on a side of the embedding rectangle, i.e. the
c   conditions 1 < IP(N) < NX and 1 < JP(N) < NY  must hold for all N.
c   Failure to ensure this will result in an error condition. In
c   addition the routine requires that there be at least one irregular
c   grid point ( NP greater than or equal to 1 ).
c
c   Boundary values at the irregular boundary points are supplied
c   through the  one-dimensional array VP, while boundary values on
c   the edge of the  embedding rectangle are supplied through the  BA,
c   BB, BC, and  BD arrays. For grid points on the edge of the embedding
c   rectangle which are not regular boundary points (i.e. not part of
c   the boundary of the solution domain), the corresponding boundary values
c   specified through the BA, BB, BC, and BD arrays are arbitrary. For
c   definiteness, they should be set to zero.
c
c   Three examples follow. In each case the grid points over
c   the rectangle are labelled according to the following
c   convention:   Interior Points - labelled 1,
c                 Regular Boundary Points - labelled 2,
c                 Irregular Boundary Points - labelled 3,
c                 Exterior Points - labelled 0.
c
c   ****  Example 1 - rectangle with a slot removed  ****
c
c    NX=17, NY=9, NP=10
c    IP=(8,8,8,8,9,10,11,11,11,11)
c    JP=(2,3,4,5,5, 5, 5, 4, 3, 2)
c
c          22222222222222222
c          21111111111111112
c          21111111111111112
c          21111111111111112
c          21111113333111112
c          21111113003111112
c          21111113003111112
c          21111113003111112
c          22222222002222222
c
c   Boundary values for the problem are specified at the irregular
c   boundary points through the VP array, and at the regular boundary
c   points through the BA, BB, BC and BD arrays. The
c   BC(I), I=1,NX array specifies boundary values on the lower edge
c   (J=1) on the embedding rectangle. These values are arbitrary for
c   I=9,10 and should be set to zero.
c
c   ****  Example 2 - rectangle with a triangular hole  ****
c
c    NX=17, NY=9, NP=20
c    IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c    JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c          22222222222222222
c          21111111111111112
c          21133333333333112
c          21113000000031112
c          21111300000311112
c          21111130003111112
c          21111113031111112
c          21111111311111112
c          22222222222222222
c
c   ****  Example 3 - triangular region  ****
c
c    NX=17, NY=9, NP=20
c    IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c    JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c          00000000000000000
c          00000000000000000
c          00033333333333000
c          00003111111130000
c          00000311111300000
c          00000031113000000
c          00000003130000000
c          00000000300000000
c          00000000000000000
c
c   There are no regular boundary points in this case, and so the
c   values of the BA, BB, BC and BD arrays are arbitrary. They
c   should be set to zero for definiteness.
c
c  Note that, in example 3, repositioning the triangular region
c  so that one of its sides coincides with one of the sides of the
c  embedding rectangle would reduce the number of irregular grid
c  points and reduce time for preprocessing and solving.
c
c  *********************************************************************
c
      REAL F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
      REAL C(NP,NP),CHARGE(NP),VP(NP),HELM,HELM2,RHO
      REAL BA(NY),BB(NY),BC(NX),BD(NX),DX,DY
      INTEGER IPS(NP),IP(NP),JP(NP)
      INTEGER IFAX(10)
      SAVE LFIRST,NX1,NY1,IFAX
      DATA LFIRST/1/

c     initialize arrays on first call or for preprocessing

      IF ((LFIRST.EQ.1).OR.(ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
        IERROR=0
        CALL CHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
        IF(IERROR.NE.0) RETURN
        NX1=NX-1
        NY1=NY-1
        CALL FAX(IFAX,NX1,4)
        CALL FFTRIG(TRIGS,NX1,4)
        LFIRST=0
      END IF

c   If required, execute the error checking routine

      IF(IERROR.EQ.-1) THEN
        IERROR=0
        CALL CHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
        IF(IERROR.NE.0) RETURN
      END IF

        HELM2=HELM*(DY**2)
        RHO=DY/DX

c     preprocessing: generate the Green's function matrix

      IF(ICFLAG.EQ.2) GOTO 5000

      IF((ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
         ICFG=ICFLAG
         CALL CPGN(C,U,WORK,TRIGS,IFAX,IPS,HELM2,
     #             RHO,NX,NY,IP,JP,NP,IERROR)
         ICFLAG=2
      ELSE
         IERROR=5
         RETURN
      END IF

      IF(ICFG.EQ.0) RETURN

5000  CONTINUE

      DO 15 J=1,NY
      DO 10 I=1,NX
         U(I,J)=F(I,J)*(DY*DY)
10    CONTINUE
15    CONTINUE

c     modify the first interior pt of the rhs to allow for
c     nonzero boundary conditions along the regular boundary

      DO 20 I=2,NX1
         U(I,2)=U(I,2)-BC(I)
         U(I,NY1)=U(I,NY1)-BD(I)
20    CONTINUE

      DO 30 J=2,NY1
         U(2,J)=U(2,J)-RHO*RHO*BA(J)
         U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
30    CONTINUE

c     get the solution in the regular domain

      CALL POISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)

c     obtain the correction to the rhs at the irregular boundary points

CVD$ NOSYNC
CDIR$ IVDEP
       DO 40 N1=1,NP
       CHARGE(N1)=VP(N1)-U(IP(N1),JP(N1))
40     CONTINUE

      CALL SGESL(C,NP,NP,IPS,CHARGE,0)

c     reconstruct the rhs

       DO 55 J=1,NY
       DO 50 I=1,NX
          U(I,J)=F(I,J)*(DY*DY)
50     CONTINUE
55     CONTINUE


      DO 60 I=2,NX1
         U(I,2)=U(I,2)-BC(I)
         U(I,NY1)=U(I,NY1)-BD(I)
60    CONTINUE

      DO 70 J=2,NY1
         U(2,J)=U(2,J)-RHO*RHO*BA(J)
         U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
70    CONTINUE

c     modify the irregular boundary points

CVD$ NOSYNC
CDIR$ IVDEP
       DO 80 N1=1,NP
          U(IP(N1),JP(N1))=U(IP(N1),JP(N1))+CHARGE(N1)
80     CONTINUE

c     get the solution in the irregular domain

      CALL POISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)

c     set the boundary conditions on the regular boundary

      DO 110 I=1,NX
         U(I,1)=BC(I)
         U(I,NY)=BD(I)
110   CONTINUE

      DO 120 J=2,NY1
         U(1,J)=BA(J)
         U(NX,J)=BB(J)
120   CONTINUE

CVD$  NOSYNC
CDIR$ IVDEP
      DO 130 N1=1,NP
         U(IP(N1),JP(N1))=VP(N1)
130   CONTINUE

      RETURN
      END

c     Routine to generate and factor the inverse capacitance matrix

      SUBROUTINE CPGN(C,R,WORK,TRIGS,IFAX,IPS,HELM,
     #                RHO,NX,NY,IP,JP,NP,IERROR)
      REAL C(NP,NP),R(NX,NY),WORK(NX,NY),TRIGS(*)
      REAL HELM,RHO
      INTEGER IP(NP),JP(NP),IPS(NP),IFAX(10)

CVD$ NOSELECT
      DO 10 N1=1,NP

      DO 25 J=1,NY
       DO 20 I=1,NX
          R(I,J)=0.
20     CONTINUE
25    CONTINUE

      I1=IP(N1)
      J1=JP(N1)
      R(I1,J1)=1.

      CALL POISS(R,HELM,RHO,WORK,TRIGS,IFAX,NX,NY)

      DO 30 N2=1,NP
         C(N2,N1)=R(IP(N2),JP(N2))
30    CONTINUE

10    CONTINUE

c     Linpack routine to factor the matrix

      CALL SGEFA(C,NP,NP,IPS,IER)

      IF(IER.NE.0) THEN
         IERROR=4
      END IF

      RETURN
      END
c
c     Routine to check input variables
c
      SUBROUTINE CHKV (NX,NY,IP,JP,NP,DX,DY,IERROR)
      INTEGER IP(NP), JP(NP)
      DATA MAXCK/40/

      NX1=NX-1
      NY1=NY-1

c     Check if the values of NX and NY are too small

      IF (NX1.LT.8.OR.NY1.LT.8) THEN
         IERROR=1
         RETURN
      END IF

c     Check if the value of NX is admissible

      LR=0
      LP=0
      LQ=0

      N1=NX-1

      DO 100 IT=1,MAXCK
         MX=MOD(N1,2)
         IF(MX.NE.0) GOTO 110
         LP=LP+1
         N1=N1/2
         IF(N1.EQ.1) GOTO 401
100   CONTINUE

110      CONTINUE

      IF(LP.EQ.0) THEN
        IERROR=2
        RETURN
      END IF

      DO 200 IT=1,MAXCK
         MX=MOD(N1,3)
         IF(MX.NE.0) GOTO 220
         LQ=LQ+1
         N1=N1/3
         IF(N1.EQ.1) GOTO 401
200   CONTINUE

220    CONTINUE

      DO 300 IT=1,MAXCK
         MX=MOD(N1,5)
         IF(MX.NE.0) GOTO 401
         LR=LR+1
         N1=N1/5
         IF(N1.EQ.1) GOTO 401
300   CONTINUE

401    CONTINUE

      NXX=(2**LP)*(3**LQ)*(5**LR) + 1

      IF(NXX.NE.NX) THEN
         IERROR=2
         RETURN
      END IF

c    Check values of the irregular boundary points

      DO 500 IN=1,NP
         IF(IP(IN).LT.2.OR.IP(IN).GT.NX-1) IERROR=3
         IF(JP(IN).LT.2.OR.JP(IN).GT.NY-1) IERROR=3
500   CONTINUE

c    Check values of DX and DY

      IF((DX.LE.0.).OR.(DY.LE.0.)) IERROR=6

      RETURN
      END
c***DECK : dcapc.f
      SUBROUTINE DCAPC(F,U,WORK,TRIGS,NX,NY,HELM,DX,DY,
     #                 C,CHARGE,IPS,VP,IP,JP,NP,ICFLAG,
     #                 BA,BB,BC,BD,IERROR)
c
c     Programmer: P. Cummins
c     Version 1.0
c
c  *********************************************************************
c
c     DCAPC is a double precision routine to solve the second-order
c     accurate finite difference approximation to the Helmholtz equation
c
c     (d/dX)(dU/dX) + (d/dY)(dU/dY) - H*U = F(X,Y)
c
c     with Dirichlet boundary conditions in irregular polygonal
c     two-dimensional domains in Cartesian coordinates using the
c     capacitance matrix method.
c
c     ******************* ON ENTRY ********************
c
c         NX   integer constant
c           The number of mesh points in the X direction
c           Note that the value of NX is restricted such
c           that (NX-1) is of the form (2**I * 3**J * 5**K),
c           where I is an integer greater than or equal to
c           one, and J and K are integers greater than or
c           equal to zero. In addition, NX-1 must be greater
c           than or equal to 8.
c
c         NY   integer constant
c           The number of mesh points in the Y direction. NY-1
c           must be greater than or equal to 8.
c
c         NP   integer constant
c           The number of irregular boundary points. NP must be
c           greater than or equal to 1.
c
c          F   double precision array (NX,NY)
c            The source function on the right hand side. The value
c            of F is arbitrary at points in the rectangle which are
c            outside the irregular solution domain. This array is
c            unchanged on exit.
c
c       HELM   double precision constant
c            The coefficient H in the Helmholtz equation given above;
c            unchanged on exit.
c
c       WORK   double precision array (NX,NY)
c            A work space array
c
c     CHARGE   double precision array (NP)
c            A work space array
c
c         DX   double precision constant
c            The discretization interval in the X direction.
c            The condition DX > 0 must hold.
c
c         DY   double precision constant
c            The discretization interval in the Y direction.
c            The condition DY > 0 must hold.
c
c         BA   double precision array (NY)
c           The boundary values on the regular boundary
c           at (I=1, J=2,NY-1).
c
c         BB   double precision array (NY)
c           The boundary values on the regular boundary
c           at (I=NX, J=2,NY-1).
c
c         BC   double precision array (nx)
c           The boundary values on the regular boundary
c           at (I=1,NX, J=1).
c
c         BD   double precision array (NX)
c           The boundary values on the regular boundary
c           at (I=1,NX, J=NY).
c
c         VP   double precision array (NP)
c           The boundary values on the irregular boundary,
c           defined by the (IP,JP) coordinates.
c
c         IP   integer array (NP)
c           The coordinates in the X direction of the
c           irregular boundary points. The condition,
c           1 < IP(N) < NX, must hold for all N.
c
c         JP   integer array (NP)
c           The coordinates in the Y direction of the
c           irregular boundary points. The condition,
c           1 < JP(N) < NY, must hold for all N.
c
c            Additional notes on defining the IP and JP arrays are
c            given below.
c
c      Note: The order in which the irregular boundary grid points, defined
c            by the (IP(N),JP(N)) pairs, are specified is immaterial.
c            However, these grid points must all be distinct; there must
c            not be a repeated (IP(N),JP(N)) pair. Failure to ensure this
c            will result in the error condition IERROR=4, defined below.
c
c     ICFLAG   integer constant
c           A flag:  ICFLAG=0 - do preprocessing only and compute
c                               the matrix C and indices IPS.
c
c                    ICFLAG=1 - do preprocessing to compute the
c                               matrix C and indices IPS, and solve
c                               the Helmholtz problem.
c
c                    ICFLAG=2 - solve the Helmholtz problem only; this
c                               assumes that preprocessing has been
c                               done and that C and IPS are defined.
c
c            Any other value for ICFLAG gives IERROR=5 and no attempt
c            is made to do preprocessing or to solve. On return,
c            ICFLAG=2.
c
c
c     ******************* ON RETURN ********************
c
c          U   double precision array (NX,NY)
c           The solution in the irregular domain, with
c           U(IP(N),JP(N)) = VP(N) for N=1,NP on the irregular
c           portion of the boundary and also with
c              U(1,J)=BA(J) for J=2,NY-1
c              U(NX,J)=BB(J) for J=2,NY-1
c              U(I,1)=BC(I) for I=1,NX
c              U(I,NY)=BD(I) for I=1,NX.
c
c          C   double precision array (NP,NP)
c           The inverse capacitance matrix. This array must not
c           be altered on successive calls to the routine with
c           different source functions when ICFLAG=2.
c
c        IPS   integer array (NP)
c           An integer vector of pivot indices. This array must
c           not be alter on successive calls with different
c           source functions when ICFLAG=2.
c
c      TRIGS   double precision array (2*(NX-1))
c            An array of coefficients required by the FFT
c            routine employed by the fast solver. This array
c            is initialized on the first call to DCAPC, or
c            whenever ICFLAG=0 or 1. The array must not be
c            altered between successive calls to the routine
c            with different source functions when ICFLAG=2.
c
c     IERROR   integer constant
c           An error flag on input variables. On return IERROR
c           has the following meaning:
c
c           IERROR=0  - normal status, no error detected.
c                 =1  - value of NX-1 or NY-1 is less than 8.
c                 =2  - illegal value for NX.
c                 =3  - illegal value for IP or JP.
c                 =4  - an error occurred in factoring the inverse
c                       capacitance matrix, i.e. the matrix is singular.
c                       This condition will occur if  the user has
c                       defined duplicate irregular boundary points.
c                 =5  - illegal value for ICFLAG.
c                 =6  - illegal value for DX or DY.
c
c     Note: Error checking is normally only done on the first call
c           to the routine. The user may force the routine to do error
c           checking at any time by setting IERROR=-1 on entry.
c
c
c   *********       NOTES ON CALLING ROUTINE DCAPC      ********
c
c     In applying the capacitance matrix method to the solution
c     of the Helmholtz equation, routine DCAPC works in two
c     stages. The first stage is a preprocessing one in which
c     the arrays C and IPS are determined. The second stage is
c     the actual solution of the problem for a given source
c     function. The flag ICFLAG determines the operation of the
c     routine:
c
c          If the routine is called with ICFLAG=0, only
c          preprocessing is done, the arrays C and IPC
c          are determined, and DCAPC immediately returns to
c          the calling program without solving. On return
c          ICFLAG=2
c
c          If the routine is called with ICFLAG=1,
c          preprocessing is done, the arrays C and IPC
c          are determined, and the solution for the
c          given source function is determined. On return
c          ICFLAG=2.
c
c          If the routine is called with ICFLAG=2,
c          the solution for the given source function
c          is determined. This assumes that preprocessing
c          has been done and that the C and IPS arrays
c          are defined. On return ICFLAG=2.
c
c     If the user wishes to solve the Helmholtz equation in a
c     given domain repeatedly with different source functions,
c     then the routine could be called with ICFLAG=1 initially.
c     Afterwards this flag is set to 2 and need not be altered.
c     If the user has stored the arrays C and IPS from a previous
c     run, then these may be supplied as input data and the
c     routine may be called initially with ICFLAG=2 to obtain
c     the solution directly.
c
c     Preprocessing need only be done once for a given domain
c     geometry and Helmholtz coefficient. The user-specified
c     variables listed below must not change between successive
c     calls to DCAPC with ICFLAG=2. If the user changes any
c     one of the following variables:
c
c                      NX,
c                      NY,
c                      DX,
c                      DY,
c                      IP,
c                      JP,
c                      NP,
c                      HELM,
c
c     then the preprocessing must be redone (ICFLAG=0 or 1). The
c     arrays TRIGS, C, and IPS should also remain unchanged on
c     successive calls with ICFLAG=2.
c
c
c   *********   NOTES ON DEFINING THE SOLUTION DOMAIN   *********
c
c   The routine solves an elliptic equation in an irregular,
c   polygonal domain. The domain is polygonal because its boundary
c   consists of straight line segments parallel or diagonal to
c   the discrete mesh. This domain is embedded in a rectangle
c   dimensioned NX by NY. The user is free to position the irregular
c   domain as desired within the embedding rectangle and the sides of
c   the rectangle may be part of the boundary of the solution domain.
c
c   The boundary of the irregular domain is described by a set of
c   grid points which are classified as 'regular' boundary points
c   and 'irregular' boundary points. A regular boundary point is one
c   which lies on one of the sides of the embedding rectangle. An
c   irregular boundary point is one which is not on a side of the
c   embedding rectangle. The routine requires that the user provide
c   the coordinates of the irregular boundary points through the
c   one-dimensional arrays IP and JP, both dimensioned NP. The user
c   must ensure that the specified coordinates of the irregular boundary
c   points do not lie on a side of the embedding rectangle, i.e. the
c   conditions 1 < IP(N) < NX and 1 < JP(N) < NY  must hold for all N.
c   Failure to ensure this will result in an error condition. In
c   addition the routine requires that there be at least one irregular
c   grid point ( NP greater than or equal to 1 ).
c
c   Boundary values at the irregular boundary points are supplied
c   through the  one-dimensional array VP, while boundary values on
c   the edge of the  embedding rectangle are supplied through the  BA,
c   BB, BC, and  BD arrays. For grid points on the edge of the embedding
c   rectangle which are not regular boundary points (i.e. not part of
c   the boundary of the solution domain), the corresponding boundary values
c   specified through the BA, BB, BC, and BD arrays are arbitrary. For
c   definiteness, they should be set to zero.
c
c   Three examples follow. In each case the grid points over
c   the rectangle are labelled according to the following
c   convention:   Interior Points - labelled 1,
c                 Regular Boundary Points - labelled 2,
c                 Irregular Boundary Points - labelled 3,
c                 Exterior Points - labelled 0.
c
c   ****  Example 1 - rectangle with a slot removed  ****
c
c    NX=17, NY=9, NP=10
c    IP=(8,8,8,8,9,10,11,11,11,11)
c    JP=(2,3,4,5,5, 5, 5, 4, 3, 2)
c
c          22222222222222222
c          21111111111111112
c          21111111111111112
c          21111111111111112
c          21111113333111112
c          21111113003111112
c          21111113003111112
c          21111113003111112
c          22222222002222222
c
c   Boundary values for the problem are specified at the irregular
c   boundary points through the VP array, and at the regular boundary
c   points through the BA, BB, BC and BD arrays. The
c   BC(I), I=1,NX array specifies boundary values on the lower edge
c   (J=1) on the embedding rectangle. These values are arbitrary for
c   I=9,10 and should be set to zero.
c
c   ****  Example 2 - rectangle with a triangular hole  ****
c
c    NX=17, NY=9, NP=20
c    IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c    JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c          22222222222222222
c          21111111111111112
c          21133333333333112
c          21113000000031112
c          21111300000311112
c          21111130003111112
c          21111113031111112
c          21111111311111112
c          22222222222222222
c
c   ****  Example 3 - triangular region  ****
c
c    NX=17, NY=9, NP=20
c    IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c    JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c          00000000000000000
c          00000000000000000
c          00033333333333000
c          00003111111130000
c          00000311111300000
c          00000031113000000
c          00000003130000000
c          00000000300000000
c          00000000000000000
c
c   There are no regular boundary points in this case, and so the
c   values of the BA, BB, BC and BD arrays are arbitrary. They
c   should be set to zero for definiteness.
c
c  Note that, in example 3, repositioning the triangular region
c  so that one of its sides coincides with one of the sides of the
c  embedding rectangle would reduce the number of irregular grid
c  points and reduce time for preprocessing and solving.
c
c  *********************************************************************
c
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
      DOUBLE PRECISION C(NP,NP),CHARGE(NP),VP(NP),HELM,HELM2,RHO
      DOUBLE PRECISION BA(NY),BB(NY),BC(NX),BD(NX),DX,DY
      INTEGER IPS(NP),IP(NP),JP(NP)
      INTEGER IFAX(10)
      SAVE LFIRST,NX1,NY1,IFAX
      DATA LFIRST/1/

c     initialize arrays on first call or for preprocessing

      IF ((LFIRST.EQ.1).OR.(ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
        IERROR=0
        CALL DCHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
        IF(IERROR.NE.0) RETURN
        NX1=NX-1
        NY1=NY-1
        CALL DFAX(IFAX,NX1,4)
        CALL DFFTRG(TRIGS,NX1,4)
        LFIRST=0
      END IF

c   If required, execute the error checking routine

      IF(IERROR.EQ.-1) THEN
        IERROR=0
        CALL DCHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
        IF(IERROR.NE.0) RETURN
      END IF

        HELM2=HELM*(DY**2)
        RHO=DY/DX

c     preprocessing: generate the Green's function matrix

      IF(ICFLAG.EQ.2) GOTO 5000

      IF((ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
         ICFG=ICFLAG
         CALL DCPGN(C,U,WORK,TRIGS,IFAX,IPS,HELM2,
     #              RHO,NX,NY,IP,JP,NP,IERROR)
         ICFLAG=2
      ELSE
         IERROR=5
         RETURN
      END IF

      IF(ICFG.EQ.0) RETURN

5000  CONTINUE

      DO 15 J=1,NY
      DO 10 I=1,NX
         U(I,J)=F(I,J)*(DY*DY)
10    CONTINUE
15    CONTINUE

c     modify the first interior pt of the rhs to allow for
c     nonzero boundary conditions along the regular boundary

      DO 20 I=2,NX1
         U(I,2)=U(I,2)-BC(I)
         U(I,NY1)=U(I,NY1)-BD(I)
20    CONTINUE

      DO 30 J=2,NY1
         U(2,J)=U(2,J)-RHO*RHO*BA(J)
         U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
30    CONTINUE

c     get the solution in the regular domain

      CALL DPOISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)

c     obtain the correction to the rhs at the irregular boundary points

CVD$ NOSYNC
CDIR$ IVDEP
       DO 40 N1=1,NP
       CHARGE(N1)=VP(N1)-U(IP(N1),JP(N1))
40     CONTINUE

      CALL DGESL(C,NP,NP,IPS,CHARGE,0)

c     reconstruct the rhs

       DO 55 J=1,NY
       DO 50 I=1,NX
          U(I,J)=F(I,J)*(DY*DY)
50     CONTINUE
55     CONTINUE


      DO 60 I=2,NX1
         U(I,2)=U(I,2)-BC(I)
         U(I,NY1)=U(I,NY1)-BD(I)
60    CONTINUE

      DO 70 J=2,NY1
         U(2,J)=U(2,J)-RHO*RHO*BA(J)
         U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
70    CONTINUE

c     modify the irregular boundary points

CVD$ NOSYNC
CDIR$ IVDEP
       DO 80 N1=1,NP
          U(IP(N1),JP(N1))=U(IP(N1),JP(N1))+CHARGE(N1)
80     CONTINUE

c     get the solution in the irregular domain

      CALL DPOISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)

c     set the boundary conditions on the regular boundary

      DO 110 I=1,NX
         U(I,1)=BC(I)
         U(I,NY)=BD(I)
110   CONTINUE

      DO 120 J=2,NY1
         U(1,J)=BA(J)
         U(NX,J)=BB(J)
120   CONTINUE

CVD$  NOSYNC
CDIR$ IVDEP
      DO 130 N1=1,NP
         U(IP(N1),JP(N1))=VP(N1)
130   CONTINUE

      RETURN
      END

c     Routine to generate and factor the inverse capacitance matrix

      SUBROUTINE DCPGN(C,R,WORK,TRIGS,IFAX,IPS,HELM,
     #                 RHO,NX,NY,IP,JP,NP,IERROR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION C(NP,NP),R(NX,NY),WORK(NX,NY),TRIGS(*)
      DOUBLE PRECISION HELM,RHO
      INTEGER IP(NP),JP(NP),IPS(NP),IFAX(10)

CVD$ NOSELECT
      DO 10 N1=1,NP

      DO 25 J=1,NY
       DO 20 I=1,NX
          R(I,J)=0.D0
20     CONTINUE
25    CONTINUE

      I1=IP(N1)
      J1=JP(N1)
      R(I1,J1)=1.D0

      CALL DPOISS(R,HELM,RHO,WORK,TRIGS,IFAX,NX,NY)

      DO 30 N2=1,NP
         C(N2,N1)=R(IP(N2),JP(N2))
30    CONTINUE

10    CONTINUE

c     Linpack routine to factor the matrix

      CALL DGEFA(C,NP,NP,IPS,IER)

      IF(IER.NE.0) THEN
         IERROR=4
      END IF

      RETURN
      END
c
c     Routine to check input variables
c
      SUBROUTINE DCHKV (NX,NY,IP,JP,NP,DX,DY,IERROR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER IP(NP), JP(NP)
      DATA MAXCK/40/

      NX1=NX-1
      NY1=NY-1

c     Check if the values of NX and NY are too small

      IF (NX1.LT.8.OR.NY1.LT.8) THEN
         IERROR=1
         RETURN
      END IF

c     Check if the value of NX is admissible

      LR=0
      LP=0
      LQ=0

      N1=NX-1

      DO 100 IT=1,MAXCK
         MX=MOD(N1,2)
         IF(MX.NE.0) GOTO 110
         LP=LP+1
         N1=N1/2
         IF(N1.EQ.1) GOTO 401
100   CONTINUE

110      CONTINUE

      IF(LP.EQ.0) THEN
        IERROR=2
        RETURN
      END IF

      DO 200 IT=1,MAXCK
         MX=MOD(N1,3)
         IF(MX.NE.0) GOTO 220
         LQ=LQ+1
         N1=N1/3
         IF(N1.EQ.1) GOTO 401
200   CONTINUE

220    CONTINUE

      DO 300 IT=1,MAXCK
         MX=MOD(N1,5)
         IF(MX.NE.0) GOTO 401
         LR=LR+1
         N1=N1/5
         IF(N1.EQ.1) GOTO 401
300   CONTINUE

401    CONTINUE

      NXX=(2**LP)*(3**LQ)*(5**LR) + 1

      IF(NXX.NE.NX) THEN
         IERROR=2
         RETURN
      END IF

c    Check values of the irregular boundary points

      DO 500 IN=1,NP
         IF(IP(IN).LT.2.OR.IP(IN).GT.NX-1) IERROR=3
         IF(JP(IN).LT.2.OR.JP(IN).GT.NY-1) IERROR=3
500   CONTINUE

c    Check values of DX and DY

      IF((DX.LE.0.D0).OR.(DY.LE.0.D0)) IERROR=6

      RETURN
      END
c***DECK : reccn.f
      SUBROUTINE RECCN(F,G,U,P,WORK,TRIGS,NX,NY,DX,DY,
     #                 BA,BB,BC,BD,ITMAX,CCON,NITT,IWRITE,IERROR)
c
c     Programmer: P. Cummins
c     Version 1.0
c
c  *********************************************************************
c
c     RECCN is a routine to solve a second-order accurate finite
c     difference approximation to the self-adjoint elliptic equation
c
c     (d/dX)(G dU/dX) + (d/dY)(G dU/dY) = F(X,Y)
c
c     with Dirichlet boundary conditions in rectangular
c     two-dimensional domains in Cartesian coordinates.
c
c     Note that the diffusion function, G, is staggered with respect
c     to the solution and the source function on the grid.
c
c         |              |             |
c  --U(I-1,J+1)-------U(I,J+1)----U(I+1,J+1)--
c         |              |             |
c         |              |             |
c         |   G(I-1,J)   |    G(I,J)   |
c         |              |             |
c         |              |             |
c  --U(I-1,J)---------U(I,J)------U(I+1,J)--
c         |              |             |
c         |              |             |
c         |  G(I-1,J-1)  |   G(I,J-1)  |
c         |              |             |
c         |              |             |
c  --U(I-1,J-1)-------U(I,J-1)----U(I-1,J+1)--
c         |              |             |
c
c
c     ******************* ON ENTRY ********************
c
c         NX   integer constant
c           The number of mesh points in the X direction
c           Note that the value of NX is restricted such
c           that (NX-1) is of the form (2**I * 3**J * 5**K),
c           where I is an integer greater than or equal to
c           one, and J and K are integers greater than or
c           equal to zero. In addition, NX-1 must be greater
c           than or equal to 8.
c
c         NY   integer constant
c           The number of mesh points in the Y direction. NY-1
c           must be greater than or equal to 8.
c
c          F   real array (NX,NY)
c            The source function on the right hand side.
c            This array is unchanged on exit.
c
c          P   real array (NX,NY)
c            The initial estimate to the solution; if the user does
c            not have an initial estimate, this array should be set
c            to zero. This array is overwritten on exit.
c
c          G   real array (NX-1,NY-1)
c            The array containing values of the diffusion function.
c            Note that this array is staggered with respect to the
c            mesh on which the source function and the solution are
c            defined. To ensure a well posed problem, the condition
c            G(I,J) > 0 must hold for all points within the rectangle.
c            Failure to ensure this will cause condition IERROR=3
c            (see below).
c
c       WORK   real array (NX,NY)
c            A work space array.
c
c         DX   real constant
c            The discretization interval in the X direction.
c            The condition DX > 0 must hold.
c
c         DY   real constant
c            The discretization interval in the Y direction.
c            The condition DY > 0 must hold.
c
c         BA   real array (NY)
c           The boundary values at (I=1, J=2,NY-1).
c
c         BB   real array (NY)
c           The boundary values at (I=NX, J=2,NY-1).
c
c         BC   real array (nx)
c           The boundary values at (I=1,NX, J=1).
c
c         BD   real array (NX)
c           The boundary values at (I=1,NX, J=NY).
c
c       CCON   real constant
c           The convergence criterion. The routine will proceed
c           with the iteration to reach a solution until the two
c           following normalized conditions are satisfied
c
c           (1)    || U(k)-U(k-1) || / ||U(k)|| <= CCON
c
c           (2)    || L(U(k))-F || / ||F|| <= CCON,
c
c           or until the number of iterations has reached the maximum
c           allowable, given by parameter ITMAX. Note that ||U(k)||
c           refers to the L2 norm of the kth iterate of U and that L
c           refers to the discrete elliptic operator. If ||F||=0, then
c           convergence is determined from condition (1) alone. CCON
c           must be greater than zero or the routine returns without
c           attempting a solution and IERROR=6.
c
c           The minimum possible value of CCON is limited by
c           the accumulation of round-off error and depends
c           mainly on the machine precision and the size of
c           the grid. The following table of approximate minimum
c           values of CCON was constructed from solution of the
c           elliptic equation, in a unit square domain with the
c           source and diffusion functions given by a uniformly
c           distributed random number, on an Alliant FX/40. This table
c           is intended as a guide to the user; the actual minimum
c           attainable value of CCON depends on the particular problem
c           and on the computer used.
c
c   (NX,NY)     Single Precision (32 bit)  Double Precision (64 bit)
c
c   (17,17)           4.E-6                      2.E-14
c
c   (33,33)           3.E-5                      5.E-14
c
c   (65,65)           2.E-4                      4.E-13
c
c   (129,129)         2.E-3                      4.E-12
c
c   (257,257)         2.E-2                      4.E-11
c
c   (513,513)         -----                      2.E-10
c
c (1025,1025)         -----                      2.E-09
c
c
c           Under normal operation, the solver will iterate until
c           the convergence conditions are satisfied and then
c           return. The convergence of the routine is monitored
c           at each iteration and if a lack of convergence is
c           detected, then the routine will return with
c           IERROR=5. Specifically the ratio D(k+1)/D(k), where
c           D(k)=|| U(k)-U(k-1) ||, is monitored for cases when
c           it is greater than one. Should this occur repeatedly
c           (three times over six iterations), then lack of
c           convergence is assumed and the routine returns to
c           the calling program with IERROR=5. This feature is
c           invoked automatically, but may be disabled by
c           setting ITMAX (see below) to a negative value.
c
c           One situation which may lead to a lack of convergence
c           arises if CCON is set to a value which is too small
c           for a given problem. In such cases, the routine
c           initially converges to the solution for a number of
c           iterations and then fails to converge further with
c           additional iterations. This may be monitored in the
c           output to unit IWRITE (see below). The remedy is to
c           increase the value of CCON or to use the double precision
c           routine DRECCN.
c
c      IWRITE  integer constant
c           Logical unit number to which the routine will write out
c           the iteration number, k, and the convergence conditions,
c
c           ANORM = || U(k)-U(k-1)|| / ||U(k)||
c
c           RESIDUAL = || L(U(k))-F || / ||F||.
c
c           This is useful for examining the iteration by iteration
c           convergence of the solver. IWRITE=0 suppresses this output.
c           If ||F||=0, then RESIDUAL is set to zero and convergence
c           depends on ANORM.
c
c      ITMAX   integer constant
c           The absolute value, |ITMAX|, gives the maximum number of
c           iteration permitted. If the number of iterations reaches
c           this maximum without obtaining a solution to the required
c           accuracy, then IERROR=4 on output. The routine accepts
c           positive or negative values for ITMAX, with negative
c           values suppressing the iterative check on convergence
c           described above under parameter CCON.
c
c     ******************* ON RETURN ********************
c
c          U   real array (NX,NY)
c           The solution in the rectangular domain, with
c           boundary conditions
c              U(1,J)=BA(J) for J=2,NY-1
c              U(NX,J)=BB(J) for J=2,NY-1
c              U(I,1)=BC(I) for I=1,NX
c              U(I,NY)=BD(I) for I=1,NX.
c
c       NITT   integer constant
c           The number of iterations required to get the
c           solution within the prescribed margin of error
c           given by CCON.
c
c      TRIGS   real array (2*(NX-1))
c            An array of coefficients required by the FFT
c            routine employed by the fast solver. This array
c            is initialized on the first call to RECCN and
c            must not be altered between successive calls to
c            the routine with different source or diffusion
c            functions.
c
c     IERROR   integer constant
c           An error flag on input variables. On return IERROR
c           has the following meaning:
c
c           IERROR=0  - normal status, no error detected.
c                 =1  - value of NX-1 or NY-1 is less than 8.
c                 =2  - illegal value for NX.
c                 =3  - an illegal value of G detected.
c                 =4  - maximum number of iterations reached.
c                 =5  - nonconvergence of the iteration detected.
c                 =6  - CCON is less than or equal to zero.
c                 =7  - illegal value for DX or DY.
c
c     Note: Error checking is normally only done on the first call
c           to the routine. The user may force the routine to do
c           error checking on subsequent calls by setting
c           IERROR=-1 on entry.
c
c  *********************************************************************
c
      PARAMETER(MMAX=6,MMAX1=MMAX-1)
      REAL F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
      REAL P(NX,NY),G(NX-1,NY-1),RHO
      REAL BA(NY),BB(NY),BC(NX),BD(NX),AN1(MMAX)
      INTEGER IFAX(10)
      SAVE LFIRST,NX1,NY1,NXSAVE,NYSAVE,IFAX
      DATA LFIRST/1/

c     initialize arrays

      IF ((LFIRST.EQ.1).OR.(NXSAVE.NE.NX).OR.(NYSAVE.NE.NY)) THEN
        IERROR=0
        CALL CHKV3(G,NX,NY,DX,DY,IERROR)
        IF(IERROR.NE.0) RETURN
        NX1=NX-1
        NY1=NY-1
        CALL FAX(IFAX,NX1,4)
        CALL FFTRIG(TRIGS,NX1,4)
        NXSAVE=NX
        NYSAVE=NY
        LFIRST=0
      END IF

c   If required, execute the error checking routine

      IF(IERROR.EQ.-1) THEN
        IERROR=0
        CALL CHKV3(G,NX,NY,DX,DY,IERROR)
        IF (IERROR.NE.0) RETURN
      END IF

      IF(CCON.LE.0.) THEN
         IERROR=6
         RETURN
      END IF


      NITT=1
      DO 5 M=1,MMAX
         AN1(M)=0.
5     CONTINUE
      RHO=DY/DX

      DO 10 I=1,NX
         P(I,1)=BC(I)
         P(I,NY)=BD(I)
10    CONTINUE

      DO 20 J=2,NY1
         P(1,J)=BA(J)
         P(NX,J)=BB(J)
20    CONTINUE

1000  CONTINUE

c     form the right hand side

      DO 35 J=2,NY1
      DO 30 I=2,NX1
      U(I,J)=(4.*DY*DY*F(I,J)-RHO*RHO*(P(I+1,J)-P(I-1,J))*
     #      (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
     #    -                  (P(I,J+1)-P(I,J-1))*
     #      (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
     #   (G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
30    CONTINUE
35    CONTINUE

c     modify the first interior pt of the rhs to allow for
c     nonzero boundary conditions along the boundary

      DO 40 I=2,NX1
         U(I,2)=U(I,2)-BC(I)
         U(I,NY1)=U(I,NY1)-BD(I)
40    CONTINUE

      DO 50 J=2,NY1
         U(2,J)=U(2,J)-RHO*RHO*BA(J)
         U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
50    CONTINUE

c     get the solution in the domain

      CALL POISS(U,0.,RHO,WORK,TRIGS,IFAX,NX,NY)

c     set the boundary conditions

      DO 120 I=1,NX
         U(I,1)=BC(I)
         U(I,NY)=BD(I)
120   CONTINUE

      DO 130 J=2,NY1
         U(1,J)=BA(J)
         U(NX,J)=BB(J)
130   CONTINUE

c     check for convergence of the solution

      ADEN=0.
      ADIFF=0.
      DL2N=0.
      DENM=0.
CVD$ ASSOC (D)
      DO 145 J=2,NY1
      DO 140 I=2,NX1
         ADEN=ADEN+U(I,J)*U(I,J)
         ADIFF=ADIFF+(U(I,J)-P(I,J))*(U(I,J)-P(I,J))

      FD=0.5*((G(I,J)+G(I,J-1))*(U(I+1,J)-U(I,J))
     #  -     (G(I-1,J)+G(I-1,J-1))*(U(I,J)-U(I-1,J)))/(DX*DX)
     #  +0.5*((G(I,J)+G(I-1,J))*(U(I,J+1)-U(I,J))
     #  -     (G(I,J-1)+G(I-1,J-1))*(U(I,J)-U(I,J-1)))/(DY*DY)

         DL2N=DL2N+(FD-F(I,J))**2
         DENM=DENM+F(I,J)**2
140   CONTINUE
145   CONTINUE

      IF(ADEN.NE.0.) THEN
        ANORM=SQRT(ADIFF/ADEN)
      ELSE
        ANORM=SQRT(ADIFF)
      END IF

      IF(DENM.NE.0.) THEN
        RRESD=SQRT(DL2N/DENM)
      ELSE IF(ADEN.NE.0.) THEN
        RRESD=SQRT(DL2N/ADEN)
      ELSE
        RRESD=SQRT(DL2N)
      END IF

c   If required, print out the iteration by iteration convergence
c   of the solver.

      IF(IWRITE.NE.0) THEN
         WRITE(IWRITE,101) NITT,ANORM,RRESD
      END IF

101   FORMAT(' ','ITERATION ',I4,' ANORM= ',1PG11.5,' RESIDUAL= ',
     #       1PG11.5)

      IF((ANORM.LE.CCON).AND.(RRESD.LE.CCON)) RETURN

      IF(ITMAX.GT.0) THEN

      DO 160 M=2,MMAX
         MM=MMAX-M+2
         AN1(MM)=AN1(MM-1)
160   CONTINUE
         AN1(1)=SQRT(ADIFF)

         RTSUM1=0.
      DO 170 M=1,MMAX1
         RAT1=0.
         IF(AN1(M+1).NE.0.) THEN
            RAT1=AN1(M)/AN1(M+1)
            IF(RAT1.GT.1.) RTSUM1=RTSUM1+1.
         END IF
170   CONTINUE

      IF((RTSUM1.GT.2.99)) THEN
         IERROR=5
         RETURN
      END IF

      END IF

      IF(NITT.GE.ABS(ITMAX)) THEN
         IERROR=4
         RETURN
      END IF

      DO 600 J=1,NY
      DO 500 I=1,NX
         P(I,J)=U(I,J)
500   CONTINUE
600   CONTINUE

      NITT=NITT+1

      GO TO 1000

      END
c
c     Routine to check input variables and the G array
c
       SUBROUTINE CHKV3 (G,NX,NY,DX,DY,IERROR)
       REAL G(NX-1,NY-1),DX,DY
       DATA MAXCK/40/

       NX1=NX-1
       NY1=NY-1

c    Check if the values of NX and NY are too small

       IF (NX1.LT.8.OR.NY1.LT.8) THEN
          IERROR=1
          RETURN
       END IF

c    Check if the value of NX is admissible

        LR=0
        LP=0
        LQ=0

        N1=NX-1

        DO 100 IT=1,MAXCK
           MX=MOD(N1,2)
           IF(MX.NE.0) GOTO 110
           LP=LP+1
           N1=N1/2
           IF(N1.EQ.1) GOTO 401
100     CONTINUE

110      CONTINUE

        IF(LP.EQ.0) THEN
          IERROR=2
          RETURN
        END IF

        DO 200 IT=1,MAXCK
           MX=MOD(N1,3)
           IF(MX.NE.0) GOTO 220
           LQ=LQ+1
           N1=N1/3
           IF(N1.EQ.1) GOTO 401
200     CONTINUE

220      CONTINUE

        DO 300 IT=1,MAXCK
           MX=MOD(N1,5)
           IF(MX.NE.0) GOTO 401
           LR=LR+1
           N1=N1/5
           IF(N1.EQ.1) GOTO 401
300     CONTINUE

401     CONTINUE

        NXX=(2**LP)*(3**LQ)*(5**LR) + 1

        IF(NXX.NE.NX) THEN
           IERROR=2
           RETURN
        END IF

        IF(IERROR.NE.0) RETURN

c   Check if the diffusion function is everywhere positive over
c   the interior of the rectangular domain

       DO 750 J=1,NY1
        DO 700 I=1,NX1
           IF(G(I,J).LE.0.) IERROR=3
700     CONTINUE
750    CONTINUE

c   Check values of DX and DY

        IF((DX.LE.0.).OR.(DY.LE.0.)) IERROR=7

        RETURN
        END
c***DECK : dreccn.f
      SUBROUTINE DRECCN(F,G,U,P,WORK,TRIGS,NX,NY,DX,DY,
     #                  BA,BB,BC,BD,ITMAX,CCON,NITT,IWRITE,IERROR)
c
c     Programmer: P. Cummins
c     Version 1.0
c
c  *********************************************************************
c
c     DRECCN is a double precision routine to solve a second-order
c     accurate finite difference approximation to the self-adjoint
c     elliptic equation
c
c     (d/dX)(G dU/dX) + (d/dY)(G dU/dY) = F(X,Y)
c
c     with Dirichlet boundary conditions in rectangular
c     two-dimensional domains in Cartesian coordinates.
c
c     Note that the diffusion function, G, is staggered with respect
c     to the solution and the source function on the grid.
c
c         |              |             |
c  --U(I-1,J+1)-------U(I,J+1)----U(I+1,J+1)--
c         |              |             |
c         |              |             |
c         |   G(I-1,J)   |    G(I,J)   |
c         |              |             |
c         |              |             |
c  --U(I-1,J)---------U(I,J)------U(I+1,J)--
c         |              |             |
c         |              |             |
c         |  G(I-1,J-1)  |   G(I,J-1)  |
c         |              |             |
c         |              |             |
c  --U(I-1,J-1)-------U(I,J-1)----U(I-1,J+1)--
c         |              |             |
c
c
c     ******************* ON ENTRY ********************
c
c         NX   integer constant
c           The number of mesh points in the X direction
c           Note that the value of NX is restricted such
c           that (NX-1) is of the form (2**I * 3**J * 5**K),
c           where I is an integer greater than or equal to
c           one, and J and K are integers greater than or
c           equal to zero. In addition, NX-1 must be greater
c           than or equal to 8.
c
c         NY   integer constant
c           The number of mesh points in the Y direction. NY-1
c           must be greater than or equal to 8.
c
c          F   double precision array (NX,NY)
c            The source function on the right hand side.
c            This array is unchanged on exit.
c
c          P   double precision array (NX,NY)
c            The initial estimate to the solution; if the user does
c            not have an initial estimate, this array should be set
c            to zero. This array is overwritten on exit.
c
c          G   double precision array (NX-1,NY-1)
c            The array containing values of the diffusion function.
c            Note that this array is staggered with respect to the
c            mesh on which the source function and the solution are
c            defined. To ensure a well posed problem, the condition
c            G(I,J) > 0 must hold for all points within the rectangle.
c            Failure to ensure this will cause condition IERROR=3
c            (see below).
c
c       WORK   double precision array (NX,NY)
c            A work space array.
c
c         DX   double precision constant
c            The discretization interval in the X direction.
c            The condition DX > 0 must hold.
c
c         DY   double precision constant
c            The discretization interval in the Y direction.
c            The condition DY > 0 must hold.
c
c         BA   double precision array (NY)
c           The boundary values at (I=1, J=2,NY-1).
c
c         BB   double precision array (NY)
c           The boundary values at (I=NX, J=2,NY-1).
c
c         BC   double precision array (nx)
c           The boundary values at (I=1,NX, J=1).
c
c         BD   double precision array (NX)
c           The boundary values at (I=1,NX, J=NY).
c
c       CCON   double precision constant
c           The convergence criterion. The routine will proceed
c           with the iteration to reach a solution until the two
c           following normalized conditions are satisfied
c
c           (1)    || U(k)-U(k-1) || / ||U(k)|| <= CCON
c
c           (2)    || L(U(k))-F || / ||F|| <= CCON,
c
c           or until the number of iterations has reached the maximum
c           allowable, given by parameter ITMAX. Note that ||U(k)||
c           refers to the L2 norm of the kth iterate of U and that L
c           refers to the discrete elliptic operator. If ||F||=0, then
c           convergence is determined from condition (1) alone. CCON
c           must be greater than zero or the routine returns without
c           attempting a solution and IERROR=6.
c
c           The minimum possible value of CCON is limited by
c           the accumulation of round-off error and depends
c           mainly on the machine precision and the size of
c           the grid. The following table of approximate minimum
c           values of CCON was constructed from solution of the
c           elliptic equation, in a unit square domain with the
c           source and diffusion functions given by a uniformly
c           distributed random number, on an Alliant FX/40. This table
c           is intended as a guide to the user; the actual minimum
c           attainable value of CCON depends on the particular problem
c           and on the computer used.
c
c   (NX,NY)     Single Precision (32 bit)  Double Precision (64 bit)
c
c   (17,17)           4.E-6                      2.E-14
c
c   (33,33)           3.E-5                      5.E-14
c
c   (65,65)           2.E-4                      4.E-13
c
c   (129,129)         2.E-3                      4.E-12
c
c   (257,257)         2.E-2                      4.E-11
c
c   (513,513)         -----                      2.E-10
c
c (1025,1025)         -----                      2.E-09
c
c
c           Under normal operation, the solver will iterate until
c           the convergence conditions are satisfied and then
c           return. The convergence of the routine is monitored
c           at each iteration and if a lack of convergence is
c           detected, then the routine will return with
c           IERROR=5. Specifically the ratio D(k+1)/D(k), where
c           D(k)=|| U(k)-U(k-1) ||, is monitored for cases when
c           it is greater than one. Should this occur repeatedly
c           (three times over six iterations), then lack of
c           convergence is assumed and the routine returns to
c           the calling program with IERROR=5. This feature is
c           invoked automatically, but may be disabled by
c           setting ITMAX (see below) to a negative value.
c
c           One situation which may lead to a lack of convergence
c           arises if CCON is set to a value which is too small
c           for a given problem. In such cases, the routine
c           initially converges to the solution for a number of
c           iterations and then fails to converge further with
c           additional iterations. This may be monitored in the
c           output to unit IWRITE (see below). The remedy is to
c           increase the value of CCON or to use the double precision
c           routine DRECCN.
c
c      IWRITE  integer constant
c           Logical unit number to which the routine will write out
c           the iteration number, k, and the convergence conditions,
c
c           ANORM = || U(k)-U(k-1)|| / ||U(k)||
c
c           RESIDUAL = || L(U(k))-F || / ||F||.
c
c           This is useful for examining the iteration by iteration
c           convergence of the solver. IWRITE=0 suppresses this output.
c           If ||F||=0, then RESIDUAL is set to zero and convergence
c           depends on ANORM.
c
c      ITMAX   integer constant
c           The absolute value, |ITMAX|, gives the maximum number of
c           iteration permitted. If the number of iterations reaches
c           this maximum without obtaining a solution to the required
c           accuracy, then IERROR=4 on output. The routine accepts
c           positive or negative values for ITMAX, with negative
c           values suppressing the iterative check on convergence
c           described above under parameter CCON.
c
c     ******************* ON RETURN ********************
c
c          U   double precision array (NX,NY)
c           The solution in the rectangular domain, with
c           boundary conditions
c              U(1,J)=BA(J) for J=2,NY-1
c              U(NX,J)=BB(J) for J=2,NY-1
c              U(I,1)=BC(I) for I=1,NX
c              U(I,NY)=BD(I) for I=1,NX.
c
c       NITT   integer constant
c           The number of iterations required to get the
c           solution within the prescribed margin of error
c           given by CCON.
c
c      TRIGS   double precision array (2*(NX-1))
c            An array of coefficients required by the FFT
c            routine employed by the fast solver. This array
c            is initialized on the first call to DRECCN and
c            must not be altered between successive calls to
c            the routine with different source or diffusion
c            functions.
c
c     IERROR   integer constant
c           An error flag on input variables. On return IERROR
c           has the following meaning:
c
c           IERROR=0  - normal status, no error detected.
c                 =1  - value of NX-1 or NY-1 is less than 8.
c                 =2  - illegal value for NX.
c                 =3  - an illegal value of G detected.
c                 =4  - maximum number of iterations reached.
c                 =5  - nonconvergence of the iteration detected.
c                 =6  - CCON is less than or equal to zero.
c                 =7  - illegal value for DX or DY.
c
c     Note: Error checking is normally only done on the first call
c           to the routine. The user may force the routine to do
c           error checking on subsequent calls by setting
c           IERROR=-1 on entry.
c
c  *********************************************************************
c
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MMAX=6,MMAX1=MMAX-1)
      DOUBLE PRECISION F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
      DOUBLE PRECISION P(NX,NY),G(NX-1,NY-1),RHO
      DOUBLE PRECISION BA(NY),BB(NY),BC(NX),BD(NX),AN1(MMAX)
      INTEGER IFAX(10)
      SAVE LFIRST,NX1,NY1,NXSAVE,NYSAVE,IFAX
      DATA LFIRST/1/

c     initialize arrays

      IF ((LFIRST.EQ.1).OR.(NXSAVE.NE.NX).OR.(NYSAVE.NE.NY)) THEN
        IERROR=0
        CALL DCHKV3(G,NX,NY,DX,DY,IERROR)
        IF(IERROR.NE.0) RETURN
        NX1=NX-1
        NY1=NY-1
        CALL DFAX(IFAX,NX1,4)
        CALL DFFTRG(TRIGS,NX1,4)
        NXSAVE=NX
        NYSAVE=NY
        LFIRST=0
      END IF

c   If required, execute the error checking routine

      IF(IERROR.EQ.-1) THEN
        IERROR=0
        CALL DCHKV3(G,NX,NY,DX,DY,IERROR)
        IF (IERROR.NE.0) RETURN
      END IF

      IF(CCON.LE.0.) THEN
         IERROR=6
         RETURN
      END IF


      NITT=1
      DO 5 M=1,MMAX
         AN1(M)=0.D0
5     CONTINUE
      RHO=DY/DX

      DO 10 I=1,NX
         P(I,1)=BC(I)
         P(I,NY)=BD(I)
10    CONTINUE

      DO 20 J=2,NY1
         P(1,J)=BA(J)
         P(NX,J)=BB(J)
20    CONTINUE

1000  CONTINUE

c     form the right hand side

      DO 35 J=2,NY1
      DO 30 I=2,NX1
      U(I,J)=(4.D0*DY*DY*F(I,J)-RHO*RHO*(P(I+1,J)-P(I-1,J))*
     #      (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
     #    -                  (P(I,J+1)-P(I,J-1))*
     #      (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
     #   (G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
30    CONTINUE
35    CONTINUE

c     modify the first interior pt of the rhs to allow for
c     nonzero boundary conditions along the boundary

      DO 40 I=2,NX1
         U(I,2)=U(I,2)-BC(I)
         U(I,NY1)=U(I,NY1)-BD(I)
40    CONTINUE

      DO 50 J=2,NY1
         U(2,J)=U(2,J)-RHO*RHO*BA(J)
         U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
50    CONTINUE

c     get the solution in the domain

      CALL DPOISS(U,0.D0,RHO,WORK,TRIGS,IFAX,NX,NY)

c     set the boundary conditions

      DO 120 I=1,NX
         U(I,1)=BC(I)