      PROGRAM PRODUIT_OMP
!
!===============================================================================
!===============================================================================
! CRIHAN     10/99
! Utilisation de directives de parallelisation OpenMP (www.openmp.org)
!
! Parallelisation implicite
!
! Calcul du vecteur R = ( U*(tU) + LAMBDA*V*(tV) )  *  ( U + V )
!
!              i.e. RESULT = MAT * VECT
!
! 1.  Initialisation distribuee du vecteur par les processus
! 2.  Initialisation distribuee des lignes de la matrice par les processus
! 3.  Calcul du produit distribue par chaque processus :
!
! Les iterations, lignes, etc ... sont distribuees sur les processus
!
!===============================================================================
!===============================================================================
!
!     DECLARATIONS
!
!===============================================================================
!===============================================================================
!
! MODULE(S) NECESSAIRE(S)
!
      use modscal
!
      IMPLICIT NONE
!
! DECLARATION DES VARIABLES
!
      REAL(rp),  ALLOCATABLE,  DIMENSION(:,:)  ::  MAT
      REAL(rp),  ALLOCATABLE,  DIMENSION(:)    ::  VECT
      REAL(rp),  ALLOCATABLE,  DIMENSION(:)    ::  RESULT
!
      REAL(rp)   ::  LAMBDA
      REAL(rp)   ::  ALPHA
      REAL(rp)   ::  BETA
      REAL(rp)   ::  scal
!
      INTEGER    ::  NBTHRDS, N
      INTEGER    ::  IERR = 0
      INTEGER    ::  ICHUNK
!
!$    EXTERNAL       OMP_GET_NUM_THREADS
!$    INTEGER    ::  OMP_GET_NUM_THREADS
!
!===============================================================================
!===============================================================================
!
!     DEBUT DU PROGRAMME
!
!===============================================================================
!===============================================================================
!
!
      OPEN ( 19, FILE = 'omp.dat', STATUS = 'OLD' )
      READ (19, *) N
      CLOSE (19)
!
      ALLOCATE ( MAT(N,N), RESULT(N), VECT(N),  STAT = IERR )
      IF (ierr /= 0) THEN
         WRITE (6,'(A)') ' Probleme allocation MAT, RESULT, VECT'
         WRITE (6,'(A)') ' IERR = ', IERR
         STOP
      END IF
!
!----- Debut de la zone parallele
!
!$OMP PARALLEL DEFAULT (SHARED)
!
!$OMP MASTER
      NBTHRDS = OMP_GET_NUM_THREADS ()
      WRITE (6,1020) 'PARALLELE : NB_THREADS = ', NBTHRDS
      ICHUNK = N / NBTHRDS
      IF ( MOD( N, NBTHRDS) /= 0 ) ICHUNK = ICHUNK + 1
 1020 FORMAT ( A, I3 )
!$OMP END MASTER
!
!$OMP BARRIER
!
!----Initialisation de la matrice et du vecteur
!
      CALL INIT_MAT_VECT (MAT, VECT, N)
!
!----Calcul du produit
!
      CALL PROD_LIGNE (MAT, VECT, RESULT, N, ICHUNK)
!
!----- Fin de la zone parallele
!
!$OMP END PARALLEL
!
!
!------Affichage des resultats
!
      WRITE (6,1030) 'Resultats', RESULT (1), RESULT (2), &
                      RESULT (N-1), RESULT (N), NBTHRDS
 1030 FORMAT ( A, 4E12.4, I4 )
!
      DEALLOCATE ( MAT, VECT, RESULT )
!
!===============================================================================
!===============================================================================
!
!     FIN DU PROGRAMME
!
!===============================================================================
!===============================================================================
!
      STOP
      END PROGRAM PRODUIT_OMP
!
!
!
!
!===============================================================
!
!     SUBROUTINE INIT_MAT_VECT : initialisation matrice vecteur
!
!     MAT       (R ,S ) Matrice du produit matriciel
!     VECT      (R ,S ) Vecteur du produit matriciel
!
!===============================================================
!
      SUBROUTINE INIT_MAT_VECT (MAT, VECT, N)
!
! MODULE(S) NECESSAIRE(S)
!
      use modscal
!
      IMPLICIT NONE
!
! DECLARATION DES ARGUMENTS
!
      INTEGER, INTENT (IN)  ::  N
!
      REAL(rp),  DIMENSION(N,N), INTENT (OUT)   ::  MAT
      REAL(rp),  DIMENSION(N),   INTENT (OUT)   ::  VECT
!
! DECLARATION DES VARIABLES LOCALES
!
      INTEGER   ::  I
      INTEGER   ::  J
      INTEGER   ::  IERR = 0
!
      REAL(rp)  ::  LAMBDA = zero
      REAL(rp)  ::  scal = 5.0_rp / 12.0_rp
!
! DEBUT DE LA SUBROUTINE
!
!
!----Initialisation des vecteurs U et V
!
!$OMP MASTER
      ALLOCATE ( U (N), V (N), STAT = IERR )
      IF (ierr /= 0) THEN
         WRITE (6,'(A)') ' Probleme allocation U, V'
         WRITE (6,'(A)') ' IERR = ', IERR
         STOP
      END IF
!$OMP END MASTER
!
!$OMP BARRIER
!
!$OMP DO PRIVATE (I)
      DO I = 1, N
         U (I) = scal
         V (I) = 0.2_rp + REAL (I)
      END DO
!$OMP END DO 
!
!
!-----Calcul de LAMBDA
!
!$OMP DO REDUCTION (+:LAMBDA) PRIVATE (I)
      DO I = 1, N
         LAMBDA = LAMBDA + U (I) * V (I)
      END DO
!$OMP END DO
!
!-----Calcul de la matrice MAT =  U*(tU) + LAMDA*V*(tV)
!$OMP DO PRIVATE (I, J)
      DO J = 1, N
!
         DO I = 1, N
            MAT (I,J) = U (I) * U (J) + LAMBDA * V (I) * V (J)
         END DO
!
      END DO
!$OMP END DO NOWAIT
!
!
!------Calcul du vecteur VECT = V + U
!$OMP DO PRIVATE (I)
      DO I = 1, N
         VECT (I) = V (I) + U (I)
      END DO
!$OMP END DO
!
!$OMP MASTER
      DEALLOCATE ( U, V )
!$OMP END MASTER
!
! FIN DE LA SUBROUTINE
!
      RETURN
      END SUBROUTINE INIT_MAT_VECT
!
!
!
!
!===============================================================
!
!     SUBROUTINE PROD_LIGNE : produit matrice vecteur parallele
!
!     MAT       (R ,E ) Matrice  du produit matriciel
!     VECT      (R ,E ) Vecteur  du produit matriciel
!     RESULT    (R, S ) Resultat du produit matriciel
!
!===============================================================
!
      SUBROUTINE PROD_LIGNE (MAT, VECT, RESULT, N, ICHUNK)
!
! MODULE(S) NECESSAIRE(S)
!
      use modscal
!
      IMPLICIT NONE
!
! DECLARATION DES ARGUMENTS
!
      INTEGER, INTENT (IN)  ::  N
      INTEGER, INTENT (IN)  ::  ICHUNK
!
      REAL(rp),  DIMENSION(N,N), INTENT (IN)   ::  MAT
      REAL(rp),  DIMENSION(N),   INTENT (IN)   ::  VECT
      REAL(rp),  DIMENSION(N),   INTENT (OUT)  ::  RESULT
!
! DECLARATION DES VARIABLES LOCALES
!
      INTEGER  ::  I
      INTEGER  ::  J
      INTEGER  ::  K
!
! DEBUT DE LA SUBROUTINE
!
!------Produit RESULT = MAT * VECT
!
!$OMP DO PRIVATE (I)
      DO I = 1, N
         RESULT (I) = zero
      END DO
!$OMP END DO
!
!$OMP DO PRIVATE (I, J, K)
      DO K = 1, N, ICHUNK
!
         DO J = 1, N
!
            DO I = K, MIN( K + ICHUNK - 1, N )
               RESULT (I) = RESULT (I) + MAT (I,J) * VECT (J)
            END DO
!
         END DO
!
      END DO
!$OMP END DO
!
! FIN DE LA SUBROUTINE
!
      RETURN
      END SUBROUTINE PROD_LIGNE

