Sample Code

The following sample code illustrates the use of PCG to solve a simple 2-d Laplace's equation on an Intel iPSC/860. For further description, see the PCG Reference Manual.

      parameter (ndim=2,nsten=5,nb=1)
      parameter (nxs=90,nys=90,nxp=4,nyp=4)
      parameter (icentr=1,inorth=2,isouth=3)
      parameter (         ieast =4,iwest =5)
      include 'fcube.h'
      include '/usr/local/pcg_fort.h'
      integer   ia(3), ja(ndim,nsten+3)
      dimension a(nsten,nxs,nys)
      dimension u(nxs,nys), b(nxs,nys)
      dimension iwk(1), fwk(1)
      dimension iparm(niparm), fparm(nfparm)
      external srigr, scg
c
c---set ia: 2-d, 5 stencil points, 1 d.o.f.
      ia(1) = ndim
      ia(2) = nsten
      ia(3) = nb
c
c---set ja
c     ---initialize to zero
      do 1 iaxis = 1, ndim
      do 1 isten = 1, nsten
 1    ja(iaxis,isten) = 0
c     ---define stencil
      ja(2,inorth ) =  1
      ja(2,isouth ) = -1
      ja(1,ieast  ) =  1
      ja(1,iwest  ) = -1
c     ---specify subgrid size      
      ja(1,nsten+1) = nxs
      ja(2,nsten+1) = nys
c     ---no. procs in each direction      
      ja(1,nsten+2) = nxp
      ja(2,nsten+2) = nyp
c     ---use gray code mapping
      ja(1,nsten+3) = GRGRAY
      ja(2,nsten+3) = GRGRAY
c
c---set a interior
      do 2 iy = 1, nys
      do 2 ix = 1, nxs
      do 3 isten = 1, nsten
 3    a(isten ,ix,iy) = -1e0
 2    a(icentr,ix,iy) =  4e0
c
c---adjust a boundaries
c     ---extract subgrid no. from
c     --- the PCG routine igrid
      ipx = mod(igrid(mynode(),
     &       ndim,ja(1,nsten+2)),nxp)
      ipy =     igrid(mynode(),
     &       ndim,ja(1,nsten+2))/nxp
c     ---zero the boundaries
      if (ipx .eq. 0  ) then
        do 11 iy = 1, nys
 11     a(iwest ,1,iy) = 0e0
      endif
      if (ipx .eq. nxp-1) then
        do 12 iy = 1, nys
 12     a(ieast ,nxs,iy) = 0e0
      endif
      if (ipy .eq. 0  ) then
        do 13 ix = 1, nxs
 13     a(isouth,ix,1) = 0e0
      endif
      if (ipy .eq. nyp-1) then
        do 14 ix = 1, nxs
 14     a(inorth,ix,nys) = 0e0
      endif
c
c---initialize i/fparm
      call sdfalt (iparm,fparm)
      iparm(LEVOUT) = LEVIT
      iparm(ITSMAX) = 1000
      fparm(ZETA)   = .001e0
c
c---set vector b
      do 20 ix = 1, nxs
      do 20 iy = 1, nys
 20   b(ix,iy) = 1e0
c
c---solve
      call spcg ( JIRT, srigr, scg,
     &  ia, ja, a, u, u, b,
     &  iwk, fwk, iparm, fparm, ier )
c
      end