PROGRAM PEAK
C**********************************************************************
C**                                                                 ***
C**   Before you install VECFEM on your computer you have to adapt  ***
C**   the include data sets to your machine. This program creates   ***
C**   the include file archi.h.                                      ***
C**                                                                 ***
C**   The program uses the include file norec.h which has to be     ***
C**   adapted before the run. Additionally the routine VEMSCD which ***
C**   is also used by VECFEM has to be adapted to your FORTRAN      ***
C**   compiler. The value DT in this program specify the accuracy   ***
C**   of the SECOND routine in seconds (e.g. DT=1.). The value      ***
C**   NMAX has to be selected to a vector length where the          ***
C**   peak performance of the general triad is nearly reached       ***
C**   (e.g. NMAX=10 000). It could happen that you have to repeat   ***
C**   the measurement with a new NMAX.                              ***
C**                                                                 ***
C**********************************************************************
C**                                                                 ***
      IMPLICIT NONE
C**                                                                 ***
C**-----------------------------------------------------------------***
C**                                                                 ***
C**      PARAMETERS :                                               ***
C**      ------------                                               ***
C**                                                                 ***
C**                    >                                            ***
      DOUBLE PRECISION  DT
      INTEGER           NMAX

      PARAMETER  (NMAX=70 000)
      PARAMETER  (DT=1.)
C**                                                                 ***
C**-----------------------------------------------------------------***
C**                                                                 ***
C**      VARIABLES AND ARRAYS :                                     ***
C**      ----------------------                                     ***
C**                                                                 ***
      DOUBLE PRECISION  X(NMAX),Y(NMAX),Z(NMAX)
      INTEGER           INDEX1(NMAX),INDEX2(NMAX)
C**
      EXTERNAL          DIA,PDIA,ROW,COL,SKY,PROW,PCOL,FROW,FCOL
C**
      INTEGER           I,ERR
      CHARACTER*10      NLOOP

      OPEN(99,FILE='peaks.h',STATUS='UNKNOWN')
C
C**** SET INDEX VECTOR :
C     ------------------
C
      DO 10 I=1,2*NMAX
        X(I)=0
        Y(I)=0
10      Z(I)=0
C
C**** START :
C     -------
C
      NLOOP='DIA '
      CALL MEASUR(NMAX,INDEX2,INDEX1,X,Y,Z,DT,DIA,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='PDIA'
      CALL MEASUR(NMAX,INDEX2,INDEX1,X,Y,Z,DT,PDIA,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='SKY '
      CALL MEASUR(NMAX,INDEX2,INDEX1,X,Y,Z,DT,SKY,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='ROW '
      CALL MEASUR(NMAX,INDEX2,INDEX1,X,Y,Z,DT,ROW,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='COL '
      CALL MEASUR(NMAX,INDEX1,INDEX2,X,Y,Z,DT,COL,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='PCOL '
      CALL MEASUR(NMAX,INDEX1,INDEX2,X,Y,Z,DT,PCOL,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='PROW '
      CALL MEASUR(NMAX,INDEX1,INDEX2,X,Y,Z,DT,PROW,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='FCOL '
      CALL MEASUR(NMAX,INDEX1,INDEX2,X,Y,Z,DT,FCOL,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999

      NLOOP='FROW '
      CALL MEASUR(NMAX,INDEX1,INDEX2,X,Y,Z,DT,FROW,NLOOP,ERR)
      IF (ERR.NE.0) GOTO 9999
9999  CONTINUE
      E    N    D
C:::::      ,,,,,MEASUR...
C
C
C
C**********************************************************************
C        1         2         3         4         5         6         7*
C**********************************************************************
C**                                                                 ***
C**                                                                 ***
      SUBROUTINE MEASUR(NMAX,INDEX1,INDEX2,X,Y,Z,DT,LOOP,NLOOP,ERR)
C**                                                                 ***
C**                                                                 ***
C**********************************************************************
C**                                                                 ***
C**      MEASUR    MEASUREMENT MVM BASE LOOP                        ***
C**                                                                 ***
C**********************************************************************
C**                                                                 ***
      IMPLICIT NONE
C**                                                                 ***
C**-----------------------------------------------------------------***
C**                                                                 ***
C**                                                                 ***
C**                    >                                            ***
      INTEGER           NMAX,ERR
      DOUBLE PRECISION  X(NMAX),Y(NMAX),Z(NMAX),DT
      INTEGER           INDEX1(NMAX),INDEX2(NMAX)
      EXTERNAL          LOOP
      CHARACTER*10      NLOOP
C**                                                                 ***
      DOUBLE PRECISION  RPEAK,R2,SIGMA,NH,NHS
      DOUBLE PRECISION  TIME1,CPU
      INTEGER           REPEAT,N,K,N0,CHECK,NHH,REP2
C
C**** START OF CALCULATION :
C     ---------------------
C
      ERR=0
C
C**** first the number of repeats is estimated :
C     ----------------------------------------
C
      REPEAT=5
      PRINT*,' '
2000  TIME1=CPU(REPEAT,1,NMAX,INDEX1,INDEX2,X,Y,Z,LOOP)
      PRINT*,'test : repeat =',REPEAT,'  time =',TIME1*REPEAT
      IF (TIME1*REPEAT.LT.DT) THEN
	 REPEAT=REPEAT*1.5
	 GOTO 2000
      ENDIF
      REPEAT=DT/TIME1
      PRINT*,'selected repeat =',REPEAT
C
C**** now the peak performance is computed :
C     -------------------------------------
C
      PRINT*,' '
      PRINT*,'peak performance :',NLOOP
      N=NMAX/2
      TIME1=CPU(REPEAT,2,N,INDEX1,INDEX2,X,Y,Z,LOOP)
      RPEAK=FLOAT(N)/TIME1*2
      PRINT*,' 1. vector length :',N
      PRINT*,'    performance [Mflops]:',RPEAK*1.D-6
C
C**** peak performance o.k. ?
C     -------------------------------------
C
      N=NMAX/4
      TIME1=CPU(REPEAT,4,N,INDEX1,INDEX2,X,Y,Z,LOOP)
      R2=FLOAT(N)/TIME1*2
      PRINT*,' 2. vector length :',N
      PRINT*,'    performance [Mflops]:',R2*1.d-6
      IF (ABS(RPEAK-R2).GT.RPEAK*.1) THEN
         WRITE(6,*) '>>increase NMAX and start again !'
         WRITE(6,*) '  abend !'
         ERR=1
         RETURN
      ENDIF
      PRINT*,'accepted !'
      RPEAK=(RPEAK+R2)/2
C
C**** n 1/2 :
C     -------------------------------------
C
      PRINT*,'n 1/2 :'
      CHECK=10
      NH=0
      NHS=0
      DO 9999 N0=11,10+CHECK
        N=N0
        K=NMAX/N
	rep2=max((1*nmax*repeat)/(K*N),1)
        TIME1=CPU(REP2,K,N,INDEX1,INDEX2,X,Y,Z,LOOP)
        NHH=TIME1*RPEAK/2-N+.5
	NH=NH+NHH
	NHS=NHS+NHH**2
	PRINT*,'    vl =',N,'   n 1/2 = ',NHH
9999  CONTINUE
      NH=NH/FLOAT(CHECK)
      SIGMA=SQRT((NHS-CHECK*NH**2)/(CHECK-1))/NH*100.
      PRINT*,'mean value ',INT(NH+.5)
      PRINT*,'varianz %',SIGMA
C
      WRITE(99,1040) NLOOP,NLOOP,NLOOP,RPEAK*1.D-6,NLOOP,INT(NH+.5)
C
C**** END OF CALCULATION
C     ------------------
C
1040  FORMAT( '       DOUBLE PRECISION   R',A4,',NH',A4
     &    /'       PARAMETER (R',A4,'=',G15.4,',NH',A4,'=',I10,')')
      R E T U R N
C-----END OF SECOND----------------------------------------------------
      E    N    D

      DOUBLE PRECISION FUNCTION CPU(REPEAT,K,N,INDEX1,
     &                              INDEX2,X,Y,Z,LOOP)
C**                                                                 ***
C**-----------------------------------------------------------------***
C**                                                                 ***
      IMPLICIT NONE
C**                                                                 ***
C**-----------------------------------------------------------------***
C**                                                                 ***
C**                    >                                            ***
      INTEGER           N,K,REPEAT
      DOUBLE PRECISION  X(N,K),Y(N,K),Z(N,K)
      INTEGER           INDEX1(N,K),INDEX2(N,K)
      EXTERNAL          LOOP
C**                                                                 ***
      DOUBLE PRECISION  TIME,VEMSCD
      INTEGER           I,J

      DO 32 I=1,K
        DO 31 J=1,N
	  INDEX1(J,I)=J
	  INDEX2(J,I)=J
31      CONTINUE
32    CONTINUE

      TIME=VEMSCD()
      DO 30 J=1,REPEAT
        DO 30 I=1,K
30        CALL LOOP(N,INDEX1(1,I),INDEX2(1,I),X(1,I),Y(1,I),Z(1,I))
      CPU=(VEMSCD()-TIME)/FLOAT(REPEAT*K)
      R E T U R N
C-----END OF CPU----------------------------------------------------
      E    N    D

      SUBROUTINE DIA(N,INDEX1,INDEX2,X,Y,Z)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      DO 10 I=1,N
10     X(I)=X(I)+Y(I)*Z(I)
      R E T U R N
      E    N    D

      SUBROUTINE PDIA(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      include"norec.h"
      DO 10 I=1,N
10     X(INDEX1(I))=X(INDEX1(I))+Y(I)*Z(INDEX1(I))
      R E T U R N
      E    N    D

      SUBROUTINE ROW(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      DO 10 I=1,N
10     X(I)=X(I)+Y(I)*Z(INDEX1(I))
      R E T U R N
      E    N    D

      SUBROUTINE COL(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      include"norec.h"
      DO 10 I=1,N
10     X(INDEX2(I))=X(INDEX2(I))+Y(I)*Z(I)
      R E T U R N
      E    N    D

      SUBROUTINE PROW(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      include"norec.h"
      DO 10 I=1,N
10     X(1)=X(1)+Y(I)*Z(INDEX2(I))
      R E T U R N
      E    N    D

      SUBROUTINE FROW(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      include"norec.h"
      DO 10 I=1,N
10     X(1)=X(1)+Y(I)*Z(I)
      R E T U R N
      E    N    D

      SUBROUTINE FCOL(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      include"norec.h"
      DO 10 I=1,N
10     X(I)=X(I)+Y(I)*Z(1)
      R E T U R N
      E    N    D

      SUBROUTINE PCOL(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      include"norec.h"
      DO 10 I=1,N
10     X(INDEX1(I))=X(INDEX1(I))+Y(I)*Z(1)
      R E T U R N
      E    N    D

      SUBROUTINE SKY(N,INDEX1,INDEX2,X,Y,Z)
      DOUBLE PRECISION    X(N),Y(N),Z(N)
      INTEGER           INDEX1(N),INDEX2(N),N,I
      include"norec.h"
      DO 10 I=1,N
10     X(INDEX2(I))=X(INDEX2(I))+Y(I)*Z(INDEX1(I))
      R E T U R N
      E    N    D