      PROGRAM NODE09
C
C     PROBLEM 9
C
C  REFERENCE:  PROBLEMS TO TEST PARALLEL AND VECTOR LANGUAGES
C              CSD-TR 516, COMPUTER SCIENCE, PURDUE UNIVERSITY
C              JOHN R. RICE, MAY 1, 1985
C
C              REVISED BY JOHN R. RICE AND J. JING, OCT. 1, 1990
C
C
C      *************************************************
C      *      Adapted for FORTRAN D benchmarking       *
C      *    by  T. HAUPT  (haupt@sccs.npac.syr.edu)    *
C      *                                               *
C      *    Northeast Parallel Architectures Center    *
C      *   at Syracuse University, Syracuse, NY, USA   *
C      *************************************************
C
C
C       VERSION MIMD/iPSC860-1.00-NODE   (F77+mp)
C       ==================================================
C
      INCLUDE 'picl_def.h'
      DOUBLE PRECISION HOST_SEC, NODE_SEC, STARTTIME
      INTEGER*4 TIMESIZE, TIMETYPE
      DATA TIMESIZE,TIMETYPE /8,10/
      INTEGER*4 ME,HOST,ALLNODES
      INTEGER*4 NPROCS,TOP,ORD,DIR
      INTEGER*4 INITTYPE
      DATA INITTYPE /1/
      DATA TOP,ORD,DIR,ALLNODES /1,1,1,-1/
C -----------------------------------------------------------
C           parameters of the initial broadcast
C -----------------------------------------------------------
      INTEGER*4 MSG(2),MSGSIZE
      DATA MSGSIZE /8/
      EQUIVALENCE (MSG(1),NK),(MSG(2),MK)
C -----------------------------------------------------------
C      parameters of messages and processes
C -----------------------------------------------------------
      INTEGER*4 PARTTYPE,ITERTYPE,ERRTYP1,ERRTYP2
      DATA PARTTYPE,ITERTYPE,ERRTYP1,ERRTYP2 /19,20,16,17/
      integer Wedge,Eedge,Sedge,Nedge
      integer WNcor,WScor,ENcor,EScor
      integer*4 mex,mey
      integer*4 nprocx,nprocy,nx,ny
      integer*4 msgW,msgE,msgN,msgS,msgWN,msgWS,msgEN,msgES
      integer*4 bytes
      data msgW,msgS,msgE,msgN /21,22,23,24/
      data msgWN,msgWS,msgEN,msgES /25,26,27,28/
      data bytes /4/

      INTEGER NK,MK
      INTEGER NINT
      REAL TSUMT,TOLER,ERR
      DATA TOLER /0.05/
      INTEGER NDIM,MDIM
      PARAMETER (NDIM=512,MDIM=512)
      INTEGER I1,J1
      REAL*4 XBUFFER(NDIM),YBUFFER(MDIM),CBUF
      REAL  U(NDIM,MDIM), T(NDIM,MDIM)
      CALL open0(NPROCS,ME,HOST)
      CALL setarc0(NPROCS,TOP,ORD,DIR)

  100 CONTINUE
      STARTTIME=clock0()
      CALL bcast0(MSG,MSGSIZE,INITTYPE,HOST)
      IF(NK.GT.0) THEN
        IF (ME.LT.NPROCS) THEN

C -----------------------------------------------------------------
C        distribute the problem
C ----------------------------------------------------------------
        WEDGE=0
        SEDGE=0
        NEDGE=0
        EEDGE=0

        CALL DISTRIB(NPROCS,NK,MK,NX,NY,NPROCX,NPROCY)

        MEY=ME/NPROCX+1
        MEX=(ME+1)-(MEY-1)*NPROCX

        if(MEX*NX.GT.NK) NX=NK-(MEX-1)*NX
        if(MEY*NY.GT.MK) NY=MK-(MEY-1)*NY

        NX2=NX+1
        NY2=NY+1

        if(MEX.EQ.1)        NEDGE=1
        if(MEX.EQ.NPROCX)   SEDGE=1
        if(MEY.EQ.1)        WEDGE=1
        if(MEY.EQ.NPROCY)   EEDGE=1


C
C     intitialization
C
         NINT=0
         ERR=1.0E10

        DO I=1,NX
           I1=I+(MEX-1)*NX
           DO J=1,NY
              J1=J+(MEY-1)*NY
              U(I+1,J+1)=I1*(I1+1)+FLOAT(J1)/(J1+1.0)
           ENDDO
        ENDDO

C
C      update U
C
       DO WHILE (ERR.GT.TOLER)

       NINT=NINT+1

C                                                 NPROCS=6
C                           North
C                                                 NPROCX=2
C                 MEY=1     MEY=2    MEY=3        NPROCY=3
C              -----------------------------
C              |         |        |        |
C       MEX=1  |    0    |   2    |   4    |
C              |         |        |        |
C West         -----------------------------        East
C              |         |        |        |
C       MEX=2  |    1    |   3    |   5    |
C              |         |        |        |
C              -----------------------------
C
C                          South
C
C         exchange messages (East <---> West)

        IF(WEDGE.EQ.0) then
          DO I=1,NX
           XBUFFER(I)=U(I+1,2)
          ENDDO
          CALL send0(XBUFFER,NX*BYTES,msgW,ME-NPROCX)
        ENDIF

        IF(EEDGE.EQ.0) then
          CALL recv0(XBUFFER,NX*BYTES,msgW)
          DO I=1,NX
            U(I+1,NY+2)=XBUFFER(I)
          ENDDO
          DO I=1,NX
           XBUFFER(I)=U(I+1,NY+1)
          ENDDO
          CALL send0(XBUFFER,NX*BYTES,msgE,ME+NPROCX)
        ENDIF

        IF(WEDGE.EQ.0) then
           CALL recv0(XBUFFER,NX*BYTES,msgE)
           DO I=1,NX
             U(I+1,1)=XBUFFER(I)
           ENDDO
        ENDIF

C         exchange messages (North <---> South)

        IF(NEDGE.EQ.0) then

          DO I=1,NY
           YBUFFER(I)=U(2,I+1)
          ENDDO
          CALL send0(YBUFFER,NY*BYTES,msgN,ME-1)
        ENDIF

        IF(SEDGE.EQ.0) then
          CALL recv0(YBUFFER,NY*BYTES,msgN)
          DO I=1,NY
            U(NX+2,I+1)=YBUFFER(I)
          ENDDO
          DO I=1,NY
           YBUFFER(I)=U(NX+1,I+1)
          ENDDO
          CALL send0(YBUFFER,NY*BYTES,msgS,ME+1)
        ENDIF

        IF(NEDGE.EQ.0) then
           CALL recv0(YBUFFER,NY*BYTES,msgS)
           DO I=1,NY
             U(1,I+1)=YBUFFER(I)
           ENDDO
        ENDIF

C       exchange messages (WN <---> ES)

        IF(WEDGE+NEDGE.EQ.0) then
          CBUF=U(2,2)
          CALL send0(CBUF,BYTES,msgWN,ME-NPROCX-1)
        ENDIF

        IF(EEDGE+SEDGE.EQ.0) then
          CALL recv0(CBUF,BYTES,msgWN)
          U(NX+2,NY+2)=CBUF
          CBUF=U(NX+1,NY+1)
          CALL send0(CBUF,BYTES,msgES,ME+NPROCX+1)
        ENDIF

        IF(WEDGE+NEDGE.EQ.0) then
          CALL recv0(CBUF,BYTES,msgES)
          U(1,1)=CBUF
        ENDIF

C       exchange messages (WS <---> EN)

        IF(WEDGE+SEDGE.EQ.0) then
          CBUF=U(NX+1,2)
          CALL send0(CBUF,BYTES,msgWS,ME-NPROCX+1)
        ENDIF

        IF(EEDGE+NEDGE.EQ.0) then
          CALL recv0(CBUF,BYTES,msgWS)
          U(1,NY+2)=CBUF
          CBUF=U(2,NY+1)
          CALL send0(CBUF,BYTES,msgEN,ME+NPROCX-1)
        ENDIF

        IF(WEDGE+SEDGE.EQ.0) then
          CALL recv0(CBUF,BYTES,msgEN)
          U(NX+2,1)=CBUF
        ENDIF


        DO i=2+NEDGE,NX2-SEDGE
         DO j=2+WEDGE,NY2-EEDGE
          t(i,j)=(u(i,j)+u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)
     &        +u(i-1,j-1)+u(i-1,j+1)+u(i+1,j-1)+u(i+1,j+1))/9.0
         ENDDO
        ENDDO

        IF(WEDGE.EQ.1) then
         DO i=2+NEDGE,NX2-SEDGE
         t(i,2)=(u(i,2)+u(i-1,2)+u(i+1,2)+u(i,3)+u(i-1,3)+u(i+1,3))/6.0
         ENDDO
        ENDIF

        IF(EEDGE.EQ.1) then
         DO i=2+NEDGE,NX2-SEDGE
         t(i,NY2)=(u(i,NY2)+u(i-1,NY2)+u(i+1,NY2)
     &            +u(i,NY2-1)+u(i-1,NY2-1)+u(i+1,NY2-1))/6.0
         ENDDO
        ENDIF

        IF(NEDGE.EQ.1) then
         DO j=2+WEDGE,NY2-EEDGE
         t(2,j)=(u(2,j)+u(2,j-1)+u(2,j+1)
     &            +u(3,j)+u(3,j-1)+u(3,j+1))/6.0

         ENDDO
        ENDIF

        IF(SEDGE.EQ.1) then
         DO j=2+WEDGE,NY2-EEDGE
         t(NX2,j)=(u(NX2,j)+u(NX2,j-1)+u(NX2,j+1)
     &            +u(NX2-1,j)+u(NX2-1,j-1)+u(NX2-1,j+1))/6.0
         ENDDO
        ENDIF

        IF(WEDGE*NEDGE.EQ.1) then
         t(2,2)=(u(2,2)+u(2,3)+u(3,3)+u(3,2))/4.0
        ENDIF

        IF(WEDGE*SEDGE.EQ.1) then
         t(NX2,2)=(u(NX2,2)+u(NX2-1,2)+u(NX2-1,3)+u(NX2,3))/4.0
        ENDIF

        IF(EEDGE*NEDGE.EQ.1) then
         t(2,NY2)=(u(2,NY2)+u(2,NY2-1)+u(3,NY2-1)+u(3,NY2))/4.0
        ENDIF

        IF(EEDGE*SEDGE.EQ.1) then
         t(NX2,NY2)=(u(NX2,NY2)+u(NX2,NY2-1)+
     &    u(NX2-1,NY2)+u(NX2-1,NY2-1))/4.0
        ENDIF

C       update matrix U
       ERR=-1.0
       DO 60 I=2,NX+1
       DO 60 J=2,NY+1
         IF(ABS(U(I,J)).GT.0.001) THEN
           ERR=AMAX1(ERR,ABS(U(I,J)-T(I,J))/U(I,J))
         ELSE
           ERR=AMAX1(ERR,ABS(U(I,J)-T(I,J))/0.001)
         ENDIF
       U(I,J)=T(I,J)
   60  CONTINUE

C        find global maximum of ERR

       CALL gmax0(ERR,1,4,ERRTYP1,0)
       CALL bcast0(ERR,4,ERRTYP2,0)
c       if(ME.EQ.0) write(6,*) nint,'global error',err

c       ERR=0
       END DO


       TSUMT=0.0
       DO I=2,NX+1
        DO J=2,NY+1
         TSUMT=TSUMT+U(I,J)
        ENDDO
       ENDDO


      CALL gsum0(TSUMT,1,4,PARTTYPE,HOST)
      CALL gmax0(NINT,1,2,ITERTYPE,HOST)
      NODE_SEC=clock0()-STARTTIME

      IF(ME.EQ.0) CALL send0(NODE_SEC,TIMESIZE,TIMETYPE,HOST)
      GOTO 100
      ENDIF

      ELSE

      CALL close0(1)
      WRITE(6,*) 'end of the job on node', ME
      STOP

      ENDIF
      END

        subroutine distrib(n,nsize,msize,nx,ny,nxp,nyp)

        implicit none
        integer n,nsize,msize,nx,ny,nxp,nyp

        if(n.eq.1) then
          nxp=1
          nyp=1
        elseif(n.eq.2) then
          if(nsize.gt.msize) then
            nxp=2
            nyp=1
           else
            nxp=1
            nyp=2
          endif
         elseif(n.eq.4) then
            nxp=2
            nyp=2
         elseif(n.eq.8) then
          if(nsize.gt.msize) then
            nxp=4
            nyp=2
           else
            nxp=2
            nyp=4
          endif
         elseif(n.eq.16) then
            nxp=4
            nyp=4
         elseif(n.eq.32) then
          if(nsize.gt.msize) then
            nxp=8
            nyp=4
           else
            nxp=4
            nyp=8
          endif
         else
            print *, 'wrong number of processors ',n
            print *, '1 processor assumed'
            nxp=1
            nyp=1
         endif

        nx=(nsize+nxp-1)/nxp
        ny=(msize+nyp-1)/nyp

         return
         end







