program iccg3d
      implicit double precision (a-h,o-z)
c   number of gridpoints (including boundary points) in x, y, and z
      parameter(nx=102,ny=102,nz=102,ndim=100)
c  nc is the leading dimension in array a (which stores the matrix)
      parameter(nc=nx*ny*nz+61208)
      common hx(nx),hy(ny),resd(ndim),x(nc),ibc(6)
      common /matrix/ a(nc,10)
      common /method/ imet
      character*80 input
      external dc,cc,mbnds
c      print 1004
c 1004 format(/' date:        )
c      print 1005
c 1005 format(/' computer: CONVEX C240, single processor, TU-Delft)
c      print 1006
c 1006 format(/'compiler:              operating system:          )
*
      print  999
  999 format(/'input the date')
      read( 5, fmt = '(a80)', end = 340 )input
      print 1004, input
 1004 format(/' date:',a)
*
      print 998
  998 format(/'input the computer')
      read( 5, fmt = '(a80)', end = 340 )input
      print 1005,input
 1005 format(/' computer:',a)
*
      print 997
  997 format(/'input the compiler/os and optimization')
      read( 5, fmt = '(a80)', end = 340 )input
      print 1006, input
 1006 format(/'compiler/operating system:',a)
  340 continue
c
********************************************************************
*  31 December 1991,  Henk A. Van Der Vorst (Utrecht University)
*                and  Jack J. Dongarra (University of Tennessee)
*
*   The purpose of this code is to measure the performance of the
*   Conjugate Gradients method for a well-structured problem. It
*   is not our goal to compare methods, nor is it our purpose to find
*   better solution methods for this particular linear system.
*   Therefore, we have adopted the following strategy.
*   50 iteration steps are carried out with the following three
*   methods:
*     standard iccg with recurrencies in the preconditioner
*     iccg with a vectorized preconditioner by a diagonal-ordering
*           scheme in x,y-plane
*     cg with diagonal scaling
*
*   The CPU-time for these 50 iteration steps is measured, as well
*   as the Mflops rate.
*
*   Those who are familiar with the algorithms and vectorization
*   approaches described in the book by Dongarra, Duff, Sorenson and
*   Van der Vorst (in particular chapter 7) will recognize these
*   elements in the Fortran source.
*   Explicit references to chapters of that book are
*   made throughout the source text of this program.
*
*   The present version of the code contains vector and parallel
comment-
*   directives for the Convex C-2 series. Also the measurement of
*   CPU-time has been done by Convex timing routines. The adaptation
*   to other computers is straight-forward.
*
*   The linear system is generated by discretizing a 3d partial
*   differential equations by finite differences over a rectangular
*   equidistant mesh over a beam, with Dirichlet boundary conditioms.
*
*   nxp, nyp and nzp denote the number of gridpoints
*   in the respective directions.
*   nxp and nyp should be less than 203, the product of
*   nxp, nyp and nzp should not exceed 1061,208. For larger
*   values the appropriate dimensions should be specified.
*
*   For other approaches for vectorizing or parallelizing (ic)cg, see
*   Dongarra, Duff, Sorenson and Van der Vorst.
*
*******************************************************************
*
c   Due to the Dirichlet boundary conditions the order of the linear
c   system is:  (nx-2)*(ny-2)*(nz-2)  (1,000,000 in this case)
c   specification of x interval:  [a1,b1]
      a1=0.d0
      b1=1.d0
c  specification of y interval:  [c1,d1]
      c1=0.d0
      d1=1.d0
c  specification of z interval:  [e1,f1]
      e1=0.d0
      f1=1.d0
c  computation of the number of unknowns in each direction
      nx2=nx-2
      ny2=ny-2
      nz2=nz-2
      n2=nx2*ny2*nz2
c some constants for discretization routine
      nxny=nx*ny
      nxnynz=nz*nxny
c  specification of the mesh sizes
      do 10 i=1,nx
 10   hx(i)=(b1-a1)/float(nx-1)
      do 20 i=1,ny
 20   hy(i)=(d1-c1)/float(ny-1)
      hzz=(f1-e1)/(nz-1.d0)
c  specification of dirichlet boundary conditions
      ibc(1)=1
      ibc(2)=1
      ibc(3)=1
      ibc(4)=1
      ibc(5)=1
      ibc(6)=1
c   eps specifies the desired error reduction.
c   nit specifies the maximum number of iteration steps
      eps=1.d-12
      nit=50
c     print 1001
c1001 format(/'  origin of linear system: 7-point finite volume',
c    0  ' discretization of'/
c    1  '  poisson equation with dirichlet bound. cond.:'/
c    2  '      in x and y direction: u=sin(t),'
c    3  ' in z-direction: u=0.'/
c    4  '      right-hand side: 0. ')
c
      print 1002,nx2,ny2,nz2,n2
 1002 format(/' number of unknowns in x-direction:',i4,
     1       /' number of unknowns in y-direction:',i4,
     2       /' number of unknowns in z-direction:',i4,
     3       /' order of the resulting matrix    :',i8)
c
c  driver routine for creation of the linear system and its solution
c
      call epsd3d(a1,b1,nx,hx,c1,d1,ny,hy,e1,f1,nz,hzz,nxny,nxnynz,
     1   dc,dc,dc,cc,cc,ibc,mbnds,mbnds,mbnds,mbnds,x,
     2   nit,eps,ierr,resd)
      stop
      end
      subroutine epsd3d(xl,xu,nx,hx,yl,yu,ny,hy,zl,zu,nz,hzz,nxny,n,
     1   dc,ec,zc,cc,fc,ibc,mbnds,mbndw,mbndn,mbnde,
     2   x,nitc,epsc,ier,resd)
      implicit double precision (a-h,o-z)
      dimension resd(1),x(*),ibc(6),hx(1),hy(1)
      common /matrix/ a(1061208,10)
      common /method/ imet
      real mflops, nflops
      external dc,ec,zc,cc,fc,mbnds,mbndw,mbndn,mbnde
      do 1000 ik=1,3
      imet=ik
      nit=nitc
      eps=epsc
      nxc=nx
      nyc=ny
      nzc=nz
      nxnyc=nxny
      nnc=n
c
c  construction of linear system. The system has a nonzero structure
c  similar to the structure in the 2D case (as in Chapter 7.2.1). The
c  nonzero diagonals of the upper triangular part of A are given by
c  a(i,1) (main diagonal), a(i,2), a(i,3) and a(i,4). The right-hand
c  side is delivered in a(i,5). For our purposes the actual values
c  of the matrix elements do not matter so much and may be replaced
c  by other values, as long as one does not take advantage of
c  knowledge of these values.
c
      call
     1 pdedi3(nxc,nyc,nzc,nxnyc,nnc,xl,xu,yl,yu,zl,zu,a(1,1),a(1,2),
     1   a(1,3),a(1,4),a(1,5),hx,hy,hzz,ibc,dc,ec,zc,cc,fc,mbnds,
     2   mbndw,mbndn,mbnde)
c
c  construction of incomplete choleski decomposition of A (of the form
c  as in Chapter 7.3. The inverse of the diagonal D is delivered in
c  a(i,6).
c
      if (imet.le.2) call minch0(a(1,6),nxc,nxnyc,nnc)
      if (imet.eq.3) then
	  do 35 i=1,nnc
 35         a(i,6)=1.d0/a(i,1)
      endif
      if (imet.eq.1) print 40
 40   format(///,'  standard iccg ')
      if (imet.eq.2) print 80
 80   format(///,'  iccg, diagonal-ordered standard preconditioning ')
      if (imet.eq.3) print 95
 95   format(///,'  cg with diagonal scaling ')
c
c  scaling of the system by the diagonal of the inc. dec.
c  this has the effect that the inc. dec. of the scaled system
c  has a unit diagonal (see chapter 7.3.1), or that for
c  unpreconditioned
c  cg the matrix has unit diagonal
c
      do 100 i=1,nnc
	 a(i,1)=a(i,1)*a(i,6)
	 a(i,10)=sqrt(a(i,6))
 100  a(i,5)=a(i,5)*a(i,10)
      do 200 i=2,nnc
 200     a(i-1,2)=a(i-1,2)*a(i,10)*a(i-1,10)
      nxc1=nxc+1
      do 300 i=nxc1,nnc
 300     a(i-nxc,3)=a(i-nxc,3)*a(i,10)*a(i-nxc,10)
      nxnyc1=nxnyc+1
      do 350 i=nxnyc1,nnc
 350     a(i-nxnyc,4)=a(i-nxnyc,4)*a(i,10)*a(i-nxnyc,10)
      do 375 i=1,nnc
 375     x(i)=0.d0
      tijd=second(t)
      call iccg(nnc,a(1,5),x(1),nxc,nyc,nzc,nit,ier,
     1   eps,resd,a(1,7),a(1,8),a(1,9))
      tijd=second(t)-tijd
c  back scaling of the solution
      do 400 i=1,nnc
 400  x(i)=x(i)*a(i,10)
c  nflops gives the number of flops per method
      nflops=((nit+1)*35)*nnc
      if (imet.eq.3) nflops=((nit+1)*22)*nnc
      mflops=nflops*1.d-6/tijd
      print 303,nit,tijd,mflops, resd(10)
 1000 continue
 303  format('  cpu time for ',i4,' iteration steps:  ',f10.3,
     1       '      mflops: ',f8.3,/
     2       '  norm of residual(10): ',d21.13)
      return
      end
      subroutine iccg(n,b,x,nx,ny,nz,nit,ierr,eps,res,
     1    r,p,ap)
      implicit double precision (a-h,o-z)
      dimension b(*),x(*),r(*),p(*),ap(*),res(*)
      common /matrix/ a(1061208,10)
      common /method/ imet
c  iccg carries out nit (ic)cg iteration steps (or less when a
c  reduction
c  by eps in the K-norm of the residual is obtained (K: preconditioner)
c  the routine is not robust, e.g., there is no test for the positive
c  definiteness of the matrix and the preconditioner. Also the stopping
c  test might be replaced by a more appropriate one.
      nxny=nx*ny
      epe=eps*eps
c computation of initial residual
      call a3x(x,r,n,nx,ny,pap)
      do 10 i=1,n
 10      r(i)=b(i)-r(i)
c preconditioning of initial residual
      if (imet.eq.1) call prec3s(r,p,n,nxny,nx,rkr)
      if (imet.eq.2) call prec3d(r,p,n,nx,ny,nz,rkr)
      if (imet.eq.3) then
	 do 15 i=1,n
 15         p(i)=r(i)
	 rkr=ddot(n,r,1,r,1)
      endif
      rrz=rkr
      res(1)=sqrt(rkr)
      nn=nit
      do 130 kk=1,nn
c computation of Ap, for given p
	 call a3x(p,ap,n,nx,ny,pap)
	 alfa=rkr/pap
c updating of approximate solution x and residual r
	 call daxpy(n,alfa,p,1,x,1)
	 call daxpy(n,-alfa,ap,1,r,1)
c computation of preconditioned residual r, as well as (r,Kr)
	 if (imet.eq.1) call prec3s(r,ap,n,nxny,nx,rkrn)
	 if (imet.eq.2) call prec3d(r,ap,n,nx,ny,nz,rkrn)
	 if (imet.eq.3) rkrn=ddot(n,r,1,r,1)
c simplified stopping test
	 res(kk+1)=sqrt(rkrn)
	 if (rkrn.lt.epe*rrz) nit=kk
	 if (kk.eq.nit) return
	 beta=rkrn/rkr
	 rkr=rkrn
	 if (imet.eq.3) then
	    do 115 i=1,n
 115           p(i)=r(i)+beta*p(i)
	 else
	    do 120 i=1,n
 120           p(i)=ap(i)+beta*p(i)
	 endif
 130  continue
      return
      end
      subroutine a3x(x, b, n, nx, ny, pap)
      implicit double precision (a-h,o-z)
      dimension x(*), b(*)
      common /method/ imet
      common /matrix/ a(1061208,10)
*
*   matrix multiplication for 7-point difference scheme.
*   for 3D problems
*
      nxny=nx*ny
      if (imet.eq.3) goto 100
      b(1)=a(1,1)*x(1)+a(1,2)*x(2)+a(1,3)*x(nx+1)+
     1      a(1,4)*x(1+nxny)
      do 10 i=2,nx
 10      b(i)=a(i-1,2)*x(i-1)+a(i,1)*x(i)+a(i,2)*x(i+1)+
     1         a(i,3)*x(i+nx)+a(i,4)*x(i+nxny)
      do 20 i=nx+1,nxny
 20      b(i)=a(i-nx,3)*x(i-nx)+a(i-1,2)*x(i-1)+a(i,1)*x(i)+
     1      a(i,2)*x(i+1)+a(i,3)*x(i+nx)+a(i,4)*x(i+nxny)
      do 30 i=nxny+1,n-nxny
	 b(i)=a(i-nxny,4)*x(i-nxny)+a(i-nx,3)*x(i-nx)+
     1        a(i-1,2)*x(i-1)+a(i,1)*x(i)+a(i,2)*x(i+1)+
     2        a(i,3)*x(i+nx)+a(i,4)*x(i+nxny)
 30   continue
      do 40 i=n-nxny+1,n-nx
 40      b(i)=a(i-nxny,4)*x(i-nxny)+a(i-nx,3)*x(i-nx)+a(i-1,2)*x(i-1)
     1   +a(i,1)*x(i)+a(i,2)*x(i+1)+a(i,3)*x(i+nx)
      do 50 i=n-nx+1,n-1
 50      b(i)=a(i-nxny,4)*x(i-nxny)+a(i-nx,3)*x(i-nx)+a(i-1,2)*x(i-1)
     1      +a(i,1)*x(i)+a(i,2)*x(i+1)
      b(n)=a(n-nxny,4)*x(n-nxny)+a(n-nx,3)*x(n-nx)+a(n-1,2)*x(n-1)+
     1      a(n,1)*x(n)
c computation of (x,Ax)
      pap=ddot(n,b,1,x,1)
      return
c computation of Ax, for given x, when A has unit main diagonal
 100  b(1)=x(1)+a(1,2)*x(2)+a(1,3)*x(nx+1)+
     1      a(1,4)*x(1+nxny)
      do 110 i=2,nx
 110     b(i)=a(i-1,2)*x(i-1)+x(i)+a(i,2)*x(i+1)+
     1         a(i,3)*x(i+nx)+a(i,4)*x(i+nxny)
      do 120 i=nx+1,nxny
 120     b(i)=a(i-nx,3)*x(i-nx)+a(i-1,2)*x(i-1)+x(i)+
     1      a(i,2)*x(i+1)+a(i,3)*x(i+nx)+a(i,4)*x(i+nxny)
      do 130 i=nxny+1,n-nxny
	 b(i)=a(i-nxny,4)*x(i-nxny)+a(i-nx,3)*x(i-nx)+
     1        a(i-1,2)*x(i-1)+x(i)+a(i,2)*x(i+1)+
     2        a(i,3)*x(i+nx)+a(i,4)*x(i+nxny)
 130  continue
      do 140 i=n-nxny+1,n-nx
 140     b(i)=a(i-nxny,4)*x(i-nxny)+a(i-nx,3)*x(i-nx)+a(i-1,2)*x(i-1)
     1   +x(i)+a(i,2)*x(i+1)+a(i,3)*x(i+nx)
      do 150 i=n-nx+1,n-1
 150     b(i)=a(i-nxny,4)*x(i-nxny)+a(i-nx,3)*x(i-nx)+a(i-1,2)*x(i-1)
     1      +x(i)+a(i,2)*x(i+1)
      b(n)=a(n-nxny,4)*x(n-nxny)+a(n-nx,3)*x(n-nx)+a(n-1,2)*x(n-1)+
     1      x(n)
c computation of innerproduct (x,Ax)
      pap=ddot(n,b,1,x,1)
      return
      end
      subroutine minch0(d,nx,nxny,n)
      implicit double precision (a-h,o-z)
      dimension d(n)
      common /matrix/ a(1061208,10)
*  modified iccg(0), 3d-version with relaxation parameter t,
*  in d the elements of the inverse of D are delivered (Chapter 7.3)
      t=0.d0
      d(1)=1./a(1,1)
      do 10 i=2,nx
 10   d(i)=1./(a(i,1)-a(i-1,2)*(a(i-1,2)+t*a(i-1,3)+t*a(i-1,4))*d(i-1))
      nx1=nx+1
      do 20 i=nx1,nxny
      inx=i-nx
 20   d(i)=1./(a(i,1)-a(i-1,2)*(a(i-1,2)+t*a(i-1,3)+t*a(i-1,4))*d(i-1)
     1        -a(inx,3)*(t*a(inx,4)+a(inx,3)+t*a(inx,2))*d(inx))
      nxny1=nxny+1
      do 30 i=nxny1,n
      inx=i-nx
      imx=i-nxny
 30   d(i)=1./(a(i,1)-a(i-1,2)*(a(i-1,2)+t*a(i-1,3)+t*a(i-1,4))*d(i-1)
     1        -a(inx,3)*(t*a(inx,2)+a(inx,3)+t*a(inx,4))*d(inx)
     2        -a(imx,4)*(t*a(imx,2)+t*a(imx,3)+a(imx,4))*d(imx))
      return
      end
      subroutine prec3s(x,px,n,nxny,nx,rkr)
      implicit double precision (a-h,o-z)
      dimension x(n),px(n)
      common /matrix/ a(1061208,10)
*   standard preconditioning iccg(0)  .
      px(1) = x(1)
      do 10 i=2,nx
	 px(i) = x(i)-a(i-1,2)*px(i-1)
   10 continue
      nx1 = nx + 1
      do 20 i=nx1,nxny
	 px(i)=x(i)-a(i-1,2)*px(i-1)-a(i-nx,3)*px(i-nx)
   20 continue
      nxny1=nxny+1
      do 25 i=nxny1,n
	 px(i)=x(i)-a(i-1,2)*px(i-1)-a(i-nx,3)*px(i-nx)-
     1          a(i-nxny,4)*px(i-nxny)
 25   continue
*
*  lower triangular part has been solved
*
c  computation of (r,Kr), where K is the preconditioner
      rkr=ddot(n,px,1,px,1)
      ij = n - 1
      do 30 i=2,nx
	 px(ij) = px(ij) - a(ij,2)*px(ij+1)
	 ij = ij - 1
   30 continue
      do 40 i=nx1,nxny
	 px(ij)=px(ij)-(a(ij,2)*px(ij+1)+a(ij,3)*px(ij+nx))
	 ij = ij - 1
   40 continue
      do 50 i=nxny1,n
	 px(ij)=px(ij)-(a(ij,2)*px(ij+1)+a(ij,3)*px(ij+nx)+
     1           a(ij,4)*px(ij+nxny))
	 ij=ij-1
 50   continue
*
*  diagonal and uppertriangular part has been solved
*
      return
      end
      subroutine prec3d(x,px,n,nx,ny,nz,rkr)
      implicit double precision (a-h,o-z)
      dimension x(n),px(n)
      common /matrix/ a(1061208,10)
      logical left
c  preconditioning for iccg
c  vectorization by diagonalwise ordering of computations
c             in (x,y) plane (Chapter 7.3.4)
      nxny=nx*ny
c solving the lower triangular part of the preconditioner
      left=.true.
      do 10 i=1,nxny
 10      px(i)=x(i)
      call precd(px(1),nxny,a(1,2),a(1,3),nx,left)
      do 30 j=2,nz
	 jl=(j-1)*nxny+1
	 ju=j*nxny
c$dir no_recurrence
	 do 20 i=jl,ju
 20         px(i)=x(i)-a(i-nxny,4)*px(i-nxny)
	 call precd(px(jl),nxny,a(jl,2),a(jl,3),nx,left)
 30   continue
      rkr=ddot(n,px,1,px,1)
      left=.false.
c solving the upper triangular part of the preconditioning
      jl=(nz-1)*nxny+1
      call precd(px(jl),nxny,a(jl,2),a(jl,3),nx,left)
      do 50 j=nz-1,1,-1
	 jl=(j-1)*nxny+1
	 ju=j*nxny
c$dir no_recurrence
	 do 40 i=jl,ju
 40         px(i)=px(i)-a(i,4)*px(i+nxny)
	 call precd(px(jl),nxny,a(jl,2),a(jl,3),nx,left)
 50   continue
      return
      end
      subroutine precd(px,nxny,a2,a3,nx,left)
      implicit double precision (a-h,o-z)
      dimension px(*),a2(*),a3(*)
      logical left
* diagonalwise computation in preconditioning over x,y-plane
* (Chapter 7.3.4)
      ny=nxny/nx
      if (.not.left) goto 15
      do 10 ii=2,nx+ny-1
	 jjmin=max(1,ii-nx+1)
	 jjmax=min(ii,ny)
	 is=(jjmin-1)*(nx-1)+ii
c$dir no_recurrence
	 do 5 j=jjmin,jjmax
	 px(is)=px(is)-a2(is-1)*px(is-1)-
     1                     a3(is-nx)*px(is-nx)
 5       is=is+nx-1
10    continue
      return
15    continue
      do 30 ii=nx+ny-1,2,-1
	 jjmin=max(1,ii-nx+1)
	 jjmax=min(ii,ny)
	 is=(jjmin-1)*(nx-1)+ii
c$dir no_recurrence
	 do 20 j=jjmin,jjmax
	 px(is)=px(is)-a2(is)*px(is+1)-
     1                     a3(is)*px(is+nx)
 20   is=is+nx-1
 30   continue
      px(1)=px(1)-a2(1)*px(2)-a3(1)*px(nx+1)
      return
      end
      subroutine pdedi3(nx,ny,nz,nxny,n,xl,xu,yl,yu,zl,zu,diag,diagx,
     1     diagy,diagz,rhs,hx,hy,hzz,ibc,dc,ec,zc,cc,fc,mbnds,mbndw,
     2     mbndn,mbnde)
      implicit double precision (a-h,o-z)
      dimension diag(*),diagx(*),diagy(*),diagz(*),rhs(*),ibc(6),
     1     hx(*),hy(*)
      external dc,ec,cc,fc,mbnds,mbndw,mbndn,mbnde,zc
*  discretisation of 3d pde.
*  mesh-sizes are assumed to be uniform
*  zc(  ) represents the coefficient of the second
*  derivative in z.
      ip=1
      n=nx*ny*nz
      do 10 i=1,n
	 diag(i)=0.
	 diagx(i)=0.
	 diagy(i)=0.
	 diagz(i)=0.
	 rhs(i)=0.
 10   continue
      z=zl
      if (ibc(5).eq.1) z=zl+hzz
      do 1000 iz=1,nz
	 nxc=nx
	 nyc=ny
	 nxnyc=nxny
	 call pdees3(nxc,nyc,nxnyc,xl,xu,yl,yu,hx,hy,mbnds,mbndw,
     1        mbndn,mbnde,ibc,dc,ec,cc,fc,diag(ip),diagx(ip),diagy(ip),
     2        rhs(ip),z)
	 call compr(nxc,nyc,nxnyc,ibc,diag(ip),diagx(ip),diagy(ip),
     1        rhs(ip))
	 ij=ip
	 zzl=z-hzz/2.
	 zzu=z+hzz/2.
	 yy=yl
	 if (ibc(1).eq.1) yy=yl+hy(1)
	 do 100 iy=1,nyc
	    dy=hy(1)
	    if (ibc(1).eq.2.and.iy.eq.1) dy=dy/2.
	    if (ibc(3).eq.2.and.iy.eq.nyc) dy=dy/2.
	    xx=xl
	    if (ibc(2).eq.1) xx=xl+hx(1)
	    do 20 ix=1,nxc
	    if (ibc(5).eq.2.and.iz.eq.1) then
	    d1=0.
	    else
	       d1=zc(xx,yy,zzl)
	    endif
	    if (ibc(6).eq.2.and.iz.eq.nz) then
	    d2=0.
	    else
	       d2=zc(xx,yy,zzu)
	    endif
	       dx=hx(1)
	       if (ibc(2).eq.2.and.ix.eq.1) dx=hx(1)/2.
	       if (ibc(4).eq.2.and.ix.eq.nxc) dx=hx(1)/2.
	       opp=dx*dy/(hzz*hzz)
	       diag (ij)=diag(ij)+opp*(d1+d2)
	       if (iz.ne.nz) diagz(ij)=-opp*d2
	       ij=ij+1
	       xx=xx+hx(1)
 20         continue
	    yy=yy+hy(1)
 100     continue
	 ip=ip+nxc*nyc
	 z=z+hzz
 1000 continue
      nx=nxc
      ny=nyc
      nxny=nxnyc
      n=nx*ny*nz
      return
      end
      subroutine mbnds(alp,bet,gam,t,z)
      implicit double precision (a-h,o-z)
*  dirichlet boundary u=sin(t), t runs from lower to upper boundary
c  limit
	 alp=1.
	 bet=0.
	 gam=sin(t)
      return
      end
      function dc(x,y,z)
      implicit double precision (a-h,o-z)
      dc=1.d0
      return
      end
      function cc(x,y,z)
      implicit double precision (a-h,o-z)
      cc=0.d0
      return
      end
      subroutine pdees3(nx, ny, nxny, a, b, c, d, hx, hy, mbnds, mbndw,
     1 mbndn, mbnde, bc, dc, ec, cc, fc, diag, sdiag, bsdiag, rhs,
     2 z)
      implicit double precision (a-h,o-z)
      dimension hx(1), hy(1), bc(4), diag(1), sdiag(1), bsdiag(1),
     1 rhs(1)
      integer bc
      external dc,ec,cc,fc,mbnds,mbndw,mbndn,mbnde
      yij = c
      ij = 0
      do 120 j=1,ny
	 h1 = 1.0
	 h3 = 1.0
	 if (j.ne.1) h1 = hy(j-1)
	 if (j.ne.ny) h3 = hy(j)
	 xij = a
	 do 110 i=1,nx
	    ij = ij + 1
	    if (i.ne.1) go to 50
	    xij = a
	    d1 = 0.
	    d2 = 0.
	    e1 = 0.
	    e2 = 0.
	    f1 = 0.
	    f2 = 0.
	    coef1 = 0.
	    coef2 = 0.
	    h2 = 1.
   50       if (j.ne.1) go to 60
	    yij = c
	    d1 = 0.
	    d4 = 0.
	    e1 = 0.
	    e4 = 0.
	    f1 = 0.
	    f4 = 0.
	    coef1 = 0.
	    coef4 = 0.
	    h1 = 1.
   60       if (i.ne.nx) go to 70
	    e3 = 0.
	    e4 = 0.
	    d3 = 0.
	    d4 = 0.
	    f3 = 0.
	    f4 = 0.
	    coef3 = 0.
	    coef4 = 0.
	    h4 = 1.
   70       if (j.ne.ny) go to 80
	    e2 = 0.
	    e3 = 0.
	    d2 = 0.
	    d3 = 0.
	    f2 = 0.
	    f3 = 0.
	    coef2 = 0.
	    coef3 = 0.
	    h3 = 1.
   80       continue
	    if (.not.(i.ne.nx .and. j.ne.ny)) go to 90
	    h4 = hx(i)
	    e3 = ec(xij+h4/2.,yij+h3/2.,z)
	    d3 = dc(xij+h4/2.,yij+h3/2.,z)
	    f3 = fc(xij+h4/2.,yij+h3/2.,z)
	    coef3 = cc(xij+h4/2.,yij+h3/2.,z)
	    diag(ij+nx) = e3
	    sdiag(ij+nx) = d3
	    bsdiag(ij+nx) = f3
	    rhs(ij+nx) = coef3
   90       if (.not.(i.ne.nx .and. j.ne.1)) go to 100
	    h4 = hx(i)
	    e4 = diag(ij)
	    d4 = sdiag(ij)
	    f4 = bsdiag(ij)
	    coef4 = rhs(ij)
  100       continue
	  diag(ij) =((h4*e4+h2*e1)/h1+(h2*e2+h4*e3)/h3+(h1*d1+h3*d2)/
     1       h2+(h3*d3+h1*d4)/h4)*0.5 + (coef1*h1*h2+coef2*h2*h3+coef3*
     2       h3*h4+coef4*h4*h1)*0.25
	    sdiag(ij) = -(h3*d3+h1*d4)*0.5/h4
	    bsdiag(ij) = -(h2*e2+h4*e3)*0.5/h3
	    rhs(ij) = (f1*h1*h2+f2*h2*h3+f3*h3*h4+f4*h4*h1)*0.25
	    xij = xij + h4
	    h2 = h4
	    e2 = e3
	    d2 = d3
	    f2 = f3
	    coef2 = coef3
	    e1 = e4
	    d1 = d4
	    f1 = f4
	    coef1 = coef4
  110    continue
	 yij = yij + h3
  120 continue
*   s-boundary.
*
      if (bc(1).eq.1) go to 140
      xij = a
       yij = c
      h2 = 1.
      h3 = hy(1)
      e2 = 0.
      do 130 i=1,nx
	 call mbnds(alp, bet, gam, xij,z)
	 h4 = 1.
	 if (i.ne.nx) h4 = hx(i)
	 e3 = 0.
	 if (i.ne.nx) e3 = ec(xij+h4/2.,yij+h3/2.,z)
	 diag(i) = diag(i) + (h2/2.*e2+h4/2.*e3)*alp/bet
	 rhs(i) = rhs(i) + (h2/2.*e2+h4/2.*e3)*gam/bet
	 e2 = e3
	 h2 = h4
	 xij = xij + h4
  130 continue
      go to 160
  140 continue
*
*   dirichlet-condition for s-boundary.
*
      ij = nx
      xij = a
      do 150 i=1,nx
	 ij = ij + 1
	 call mbnds(alp, bet, gam, xij,z)
	 rhs(ij) = rhs(ij) - bsdiag(i)*gam
	 if (i.eq.nx) go to 150
	 xij = xij + hx(i)
  150 continue
  160 continue
*
*   w-boundary.
*
      xij = a
      yij = c
      d4 = 0.0
      h1 = 1.0
      h4 = hx(1)
      if (bc(2).eq.1) go to 180
      iw = 1
      do 170 i=1,ny
	 call mbndw(alp, bet, gam, yij,z)
	 h3 = 1.
	 if (i.ne.ny) h3 = hy(i)
	 d3 = 0.0
	 if (i.ne.ny) d3 = dc(xij+h4/2.,yij+h3/2.,z)
	 diag(iw) = diag(iw) + (h1/2.*d4+h3/2.*d3)*alp/bet
	 rhs(iw) = rhs(iw) + (h1/2.*d4+h3/2.*d3)*gam/bet
	 iw = iw + nx
	 d4 = d3
	 h1 = h3
	 yij = yij + h3
  170 continue
      go to 200
  180 continue
*
*   dirichlet-condition for w-boundary.
*
      iw = 2
      do 190 i=1,ny
	 call mbndw(alp, bet, gam, yij,z)
	 rhs(iw) = rhs(iw) - sdiag(iw-1)*gam
	 iw = iw + nx
	 if (i.eq.ny) go to 190
	 yij = yij + hy(i)
  190 continue
  200 continue
*
* the north-boundary..
*
      xij = a
      yij = d
      h1 = hy(ny-1)
      e1 = 0.0
      h2 = 1.0
      if (bc(3).eq.1) go to 220
      ij = nx*(ny-1) + 1
      do 210 i=1,nx
	 call mbndn(alp, bet, gam, xij,z)
	 h4 = 1.
	 if (i.ne.nx) h4 = hx(i)
	 e4 = 0.0
	 if (i.ne.nx) e4 = ec(xij+h4/2.,yij-h1/2.,z)
	 diag(ij) = diag(ij) + (h2/2.*e1+h4/2.*e4)*alp/bet
	 rhs(ij) = rhs(ij) + (h2/2.*e1+h4/2.*e4)*gam/bet
	 ij = ij + 1
	 h2 = h4
	 e1 = e4
	 xij = xij + h4
  210 continue
      go to 240
  220 continue
      ij = nx*(ny-2)
      do 230 i=1,nx
	 ij = ij + 1
	 call mbndn(alp, bet, gam, xij, z)
	 rhs(ij) = rhs(ij) - bsdiag(ij)*gam
	 bsdiag(ij) = 0.0
	 if (i.eq.nx) go to 230
	 xij = xij + hx(i)
  230 continue
  240 continue
*
* east-boundary..
*
      xij = b
      yij = c
      d1 = 0.0
      h1 = 1.0
      h2 = hx(nx-1)
      if (bc(4).eq.1) go to 260
      ie = nx
      do 250 i=1,ny
	 call mbnde(alp, bet, gam, yij, z)
	 h3 = 1.
	 if (i.ne.ny) h3 = hy(i)
	 d2 = 0.0
	 if ((i.ne.ny) .or. ((bc(3).eq.1) .and. (bc(1).ne.3))) d2 =
     1    dc(xij-h2/2.,yij+h3/2.,z)
	 diag(ie) = diag(ie) + (h1/2.*d1+h3/2.*d2)*alp/bet
	 rhs(ie) = rhs(ie) + (h1/2.*d1+h3/2.*d2)*gam/bet
	 ie = ie + nx
	 d1 = d2
	 h1 = h3
	 yij = yij + h3
  250 continue
      go to 290
*  dirichlet-condition on east-boundary..
*
  260 continue
*
      ie = nx - 1
      do 270 i=1,ny
	 call mbnde(alp, bet, gam, yij, z)
	 rhs(ie) = rhs(ie) - sdiag(ie)*gam
	 sdiag(ie) = 0.0
	 ie = ie + nx
	 if (i.eq.ny) go to 270
	 yij = yij + hy(i)
  270 continue
  290 continue
      nxny = nx*ny
      return
      end
      subroutine compr(nx, ny, nxny, ibc, diag, sdiag, bsdiag, rhs)
      implicit double precision (a-h,o-z)
      dimension ibc(4), diag(1), sdiag(1), bsdiag(1), rhs(1)
      if (ibc(1).ne.1) go to 20
*
*   elimination of south-boundary. system becomes nx smaller.
*
      ny = ny - 1
      nxny = nx*ny
      do 10 i=1,nxny
	 inx = i + nx
	 diag(i) = diag(inx)
	 sdiag(i) = sdiag(inx)
	 bsdiag(i) = bsdiag(inx)
	 rhs(i) = rhs(inx)
   10 continue
   20 continue
      if (ibc(2).ne.1) go to 70
*
* west-boundary eliminated, system becomes ny smaller..
*
      isch = 1
      nx = nx - 1
      ij = 1
      do 60 i=1,ny
	 do 50 j=1,nx
	    ijisch = ij + isch
	    diag(ij) = diag(ijisch)
	    sdiag(ij) = sdiag(ijisch)
	    bsdiag(ij) = bsdiag(ijisch)
	    rhs(ij) = rhs(ijisch)
	    ij = ij + 1
   50    continue
	 isch = isch + 1
   60 continue
   70 continue
      if (.not.((ibc(1).eq.3) .or. (ibc(3).eq.1))) go to 100
*
* dirichlet-condition on north-boundary..
*
      ny = ny - 1
  100 continue
      if (.not.((ibc(2).eq.3) .or. (ibc(4).eq.1))) go to 130
*
* east-boundary is eliminated, system becomes ny smaller..
*
      nx = nx - 1
      isch = 1
      ij = nx + 1
      do 120 i=2,ny
	 do 110 j=1,nx
	    ijisch = ij + isch
	    diag(ij) = diag(ijisch)
	    sdiag(ij) = sdiag(ijisch)
	    bsdiag(ij) = bsdiag(ijisch)
	    rhs(ij) = rhs(ijisch)
	    ij = ij + 1
  110    continue
	 isch = isch + 1
  120 continue
  130 continue
      nxny = nx*ny
      return
      end
      function second(t)
      implicit double precision (a-h,o-z)
c second(t) should deliver the accumulated CPU-time in seconds.
c to be replaced by appropriate function or call
      second=cputime(0.)
      return
      end





.