      PROGRAM PROB09
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 SIMD/DECmpp-1.00
C       ==================================================
C
      INTEGER KASES,NK,MK
      PARAMETER (KASES=4)
      INTEGER N(KASES),M(KASES)
      DATA N /64,256,256,512/
      DATA M /128,128,256,512/
      INTEGER ITIME
      REAL TIME


      DO 50 K = 1, KASES

      CALL mpTIMERstart()
      DO MANY=1,3
      NK=N(K)
      MK=M(K)
      CALL DOIT(NK,MK,TSUMT,NINT)
      ENDDO
      ITIME=mpTIMERelapsed()

      PRINT 60, N(K),M(K)
  60  FORMAT ('PROBLEM 9 WITH N,M =',I6,2X,I6)
      PRINT *,'GIVES TSUMT =',TSUMT
      PRINT *,'AFTER ',NINT,' ITERATIONS'

      TIME=FLOAT(ITIME)/1000
      PRINT *, 'DECmpp ELAPSED TIME:',TIME,' SECONDS'

   50 CONTINUE

      STOP
      END


      SUBROUTINE DOIT(NK,MK,TSUMT,NINT)
      INTEGER NK,MK
      INTEGER NINT
      REAL TSUMT,TOLER,ERR
      DATA TOLER /0.05/
      REAL, ARRAY(NK) :: TN
      REAL, ARRAY(MK) :: TM
      REAL, ARRAY(NK,MK) :: T,RIP,RIN,PATTERN,U,ERRM
C
C     intitialization
C
         NINT=0
         ERR=1.0E10

       PATTERN=4.0
       PATTERN(2:NK-1,2:MK-1)=9.0
       PATTERN(1,2:MK-1)=6.0
       PATTERN(NK,2:MK-1)=6.0
       PATTERN(2:NK-1,1)=6.0
       PATTERN(2:NK-1,MK)=6.0

       TN=[1:NK]*([1:NK]+1)
       TM=REAL([1:MK])/([1:MK]+1.0)
       T=SPREAD(TN,2,MK)+SPREAD(TM,1,NK)
C
C      update U
C
       DO WHILE (ERR.GT.TOLER)

       NINT=NINT+1
       U=T
       T=0

       RIP=EOSHIFT(U,DIM=1,SHIFT=1)
       RIN=EOSHIFT(U,DIM=1,SHIFT=-1)

       T=U+RIP+RIN
       T=T+EOSHIFT(U,DIM=2,SHIFT=1)
       T=T+EOSHIFT(U,DIM=2,SHIFT=-1)
       T=T+EOSHIFT(RIP,DIM=2,SHIFT=1)
       T=T+EOSHIFT(RIP,DIM=2,SHIFT=-1)
       T=T+EOSHIFT(RIN,DIM=2,SHIFT=1)
       T=T+EOSHIFT(RIN,DIM=2,SHIFT=-1)


       T=T/PATTERN

       WHERE (ABS(U).GT.0.001)
          ERRM=ABS((U-T)/U)
       ELSEWHERE
          ERRM=ABS((U-T)/0.001)
       ENDWHERE
       ERR=MAXVAL(ERRM)
c
c       U=T


       END DO


       TSUMT= SUM(U(1:NK,1:MK))

      RETURN
      END





