program main
  !  Programme test de perfo tableaux f90
  !   Reza Ansari - Avril 2003
  implicit none
  parameter NR = 1000
  parameter NC = 500
  real(8) , dimension(:,:), allocatable :: mx1, mx2, mx3, mx4, mx5
  integer i,j,k,ntime
  real(8) x,y,z,c1,c2,c3
  integer count
  integer nra,nca

  count = 0
  ntime = 50 
  nra = NR
  nca = NC

  allocate(mx1(NR,NC))
  allocate(mx2(NR,NC))
  allocate(mx3(NR,NC))
!  allocate(mx4(NR,NC))
  allocate(mx5(NR,NC))

  call finittim()
  print *,' f90mtx - NR=' , NR, ' NC=' , NC
  print *,' Remplissage MX1 MX2 (NRxNC) - NTimes=', ntime
  do k=1,ntime/5
     x = k*0.00314
     y = sin(X)
     do i=1,NR
        do j=1,NC
           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) 
     call random_number(x)
     call random_number(y)
     call random_number(z)
     i = x*NR
     j = y*NC
     if ((z>0.98).or.(k.eq.15).or.(k.eq.35)) &
          print *, 'FillMtx/Check: ', & 
          k, mx1(i,j), mx2(i,j), mx3(i,j), &
          ' x,y,z=',z,y,z
  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 opemtx(nra, nca, c1,c2,c3, mx1,mx2,mx3,mx5,count)
     call random_number(x)
     call random_number(y)
!     call random_number(z)
     z = 0
     i = x*NR
     j = y*NC
     if ((z>0.98).or.(k.eq.15).or.(k.eq.35))   &
          print *, 'FillMtx/Check: ', k, i,j, &
          mx1(i,j),mx2(i,j), mx3(i,j), &
          '  c1,c2,c3=', c1,c2,c3, & 
          '  x,y,z=' , x,y,z
     
  enddo
  print *, 'FIN  MX5 = C1*MX1 + MX1 * MX2*C2 + MX2C3 ' 
  call fprttim(2)

  print *, '   NCount appel opemtx = ' , count 
  print *, ' ------ FIN programme f90mtx ---------- ' 
end program main

subroutine opemtx(nra,nca,c1,c2,c3,mx1,mx2,mx3,mx,count)
  integer  :: nra,nca,count
  real(8), dimension(nra,nca) :: mx1,mx2,mx3
  real(8), dimension(nra,nca) :: mx
  real(8)  :: c1,c2,c3
  mx = c1*mx1+c2*mx2+c3*mx3
  count = count+1
  if ((count.eq.9).or.(count.eq.36)) &
       print *, 'opemtx-count=' , count, ' mx1,2,3= ', &
       mx1(count,19), mx2(count,19), mx3(count, 19), &
       ' c1,c2,c3=', c1,c2,c3, &
       ' mx= ' , mx(count,19)
  return
end subroutine opemtx