PROGRAM MAIN
      IMPLICIT NONE
C      PARAMETER M = 50000
      DOUBLE PRECISION MX1(1000,500), MX2(1000,500), MX3(1000,500)
      DOUBLE PRECISION MX4(1000,500), MX5(1000,500)
      INTEGER I,J,K,NTIME,TYP
      DOUBLE PRECISION X,Y,Z,C1,C2,C3

      CALL FInitTim()   
      NTIME = 50 
      TYP = 3
      PRINT *,' Remplissage MX1 MX2 (1000x500) - NTimes=', NTIME
      DO K=1,NTIME
        X = FLOAT(K)*0.00314
        Y = SIN(X)
        CALL FILLMTX(K,X,Y,MX1,MX2,MX3)
      ENDDO

      PRINT *, 'FIN Remplissage ' 
      CALL FPrtTim(1)

      DO K=1,NTIME
        C1 = K*0.035
        C2 = COS(C1)+0.04
        C3 = SIN(C1*5.)
        CALL OPERMTX(K, MX1, MX2, MX5, C1, C2, C3)
      ENDDO
      PRINT *, 'FIN  MX5 = C1*MX1 + MX1 * MX2*C2 + MX2C3 ' 
      CALL FPrtTim(2)

      DO K=1,NTIME
        C1 = K*0.055
        C2 = COS(C1*1.6)+0.023
        C3 = SIN(C1*2.5)
        CALL OPADDMTX(K, MX1, MX2, MX3, MX5, C1, C2, C3)
      ENDDO
      PRINT *, 'FIN  MX5 = C1*MX1 + C2*MX2 + C3*MX3 ' 
      CALL FPrtTim(3)

      DO K=1,NTIME
        C1 = K*0.055
        C2 = COS(C1*1.6)+0.023
        C3 = SIN(C1*2.5)
        CALL MTXBLAS(K, MX1, MX2, MX3, MX4, MX5, C1, C2, C3)
      ENDDO
      PRINT *, 'FIN  MTXBLAS MX5 = C1*MX1 + C2*MX2 + C3*MX3 ' 
      CALL FPrtTim(4)

      DO K=1,NTIME
        C1 = K*0.135
        C2 = COS(6.*C1)+0.04
        C3 = SIN(C1*8.5)
        CALL OPEFMTX(K, MX1, MX2, MX5, C1, C2, C3)
      ENDDO
      PRINT *, 'FIN  MX5 = C1*MX1 + MX3 + COS(C2*MX2*MX1) ' 
      CALL FPrtTim(5)

      PRINT *, ' ------ FIN programme ---------- ' 
      END

      SUBROUTINE FILLMTX(K, X,Y, MX1, MX2, MX3)
      DOUBLE PRECISION X,Y
      DOUBLE PRECISION MX1(1000,500), MX2(1000,500), MX3(1000,500)
      INTEGER I,J,K
      DO I=1,1000
         DO J=1,500
            MX1(I,J) = K*300.+10.*I+J+X
            MX2(I,J) = K*550.+20.*I+2.*J+Y
            MX3(I,J) = K*860.+40.*I+7.*J+Y*3.14
         ENDDO
      ENDDO
      IF (K.EQ.5) PRINT *, 'FILLMTX - 5 -> ', MX1(12,14) 
      IF (K.EQ.15) PRINT *, 'FILLMTX - 15 -> ', MX2(12,15) 

      RETURN
      END

      SUBROUTINE OPERMTX(K, MX1, MX2, MX5, C1, C2, C3)
      INTEGER I,J,K
      DOUBLE PRECISION MX1(1000,500), MX2(1000,500)
      DOUBLE PRECISION MX5(1000,500)
      DOUBLE PRECISION C1,C2,C3
      DO I=1,1000
         DO J=1,500
            MX5(I,J) = C1*MX1(I,J)+MX1(I,J)*(MX2(I,J)*C2)+MX2(I,J)*C3
         ENDDO
      ENDDO
      IF (K.EQ.5) PRINT *, 'OPERMTX - 5 -> ', MX5(20,5)
      IF (K.EQ.15) PRINT *, 'OPEEMTX - 15 -> ', MX5(20,15)
      RETURN
      END

      SUBROUTINE OPEFMTX(K, MX1, MX2, MX5, C1, C2, C3)
      INTEGER I,J,K
      DOUBLE PRECISION MX1(1000,500), MX2(1000,500)
      DOUBLE PRECISION MX5(1000,500)
      DOUBLE PRECISION C1,C2,C3
      DO I=1,1000
         DO J=1,500
            MX5(I,J) = C1*MX1(I,J)+MX2(I,J)+
     +                 C3*COS(C2*MX2(I,J)*MX1(I,J))
         ENDDO
      ENDDO
      IF (K.EQ.5) PRINT *, 'OPEFMTX/COS - 5 -> ', MX5(20,5)
      IF (K.EQ.15) PRINT *, 'OPEFMTX/COS - 15 -> ', MX5(20,15)
      RETURN
      END

      SUBROUTINE OPADDMTX(K, MX1, MX2, MX3, MX5, C1, C2, C3)
      INTEGER I,J,K
      DOUBLE PRECISION MX1(1000,500), MX2(1000,500)
      DOUBLE PRECISION MX3(1000,500), MX5(1000,500)
      DOUBLE PRECISION C1,C2,C3
      DO I=1,1000
         DO J=1,500
            MX5(I,J) = C1*MX1(I,J)+C2*MX2(I,J)+MX3(I,J)*C3
         ENDDO
      ENDDO
      IF (K.EQ.5) PRINT *, 'OPADDMTX - 5 -> ', MX5(20,5)
      IF (K.EQ.15) PRINT *, 'OPADDMTX - 15 -> ', MX5(20,15)

      RETURN
      END

      SUBROUTINE MTXBLAS(K, MX1, MX2, MX3, MX4, MX5, C1, C2, C3)
      INTEGER I,J,K,NTOT
      DOUBLE PRECISION MX1(1000,500), MX2(1000,500)
      DOUBLE PRECISION MX5(1000,500)
      DOUBLE PRECISION C1,C2,C3
      NTOT = 1000*500
C  Sur OSF :
C      CALL DVCAL(NTOT, C1, MX1, 1, MX5, 1)
C      CALL DZAXPY(NTOT, C2, MX2, 1, MX5, 1, MX4, 1)
C      CALL DZAXPY(NTOT, C3, MX3, 1, MX4, 1, MX5, 1)
C  Sinon 
      CALL DCOPY(NTOT, MX1, 1, MX5, 1)
      CALL DSCAL(NTOT, C1, MX5, 1)
      CALL DAXPY(NTOT, C2, MX2, 1, MX5, 1)
      CALL DAXPY(NTOT, C3, MX3, 1, MX3, 1)
      IF (K.EQ.5) PRINT *, 'MTXBLAS - 5 -> ', MX5(20,5),
     + ' NTot=' , NTOT
      IF (K.EQ.15) PRINT *, 'MTXBLAS - 15 -> ', MX5(20,15), 
     + ' NTot=' , NTOT
      RETURN
      END