!************************************************************************!
!Resolution de l'equation de Helmholtz en utilisant le Gradient Conjugue
!     sur le domaine [0,1]x[0,1]x[0,1] par une methode aux differences finies,
!     avec comme solveur le Gradient Conjugue.
!     u - \lambda * Delta  u  =  f(x,y,z)
!                             =  2 * { z*(z-1) * [ x*(x-1) + y*(y-1) ]
!                                    + x*(x-1)*y*(y-1) }
!     u sur les bords vaut 0
!     La solution exacte est u =  x*y*z*(x-1)*(y-1)*(z-1)
!
!     Dans cette version, on se donne
!     le nombre de points interieurs total suivant x (ntx),
!     le nombre de points inferieurs total suivant y (nty) et 
!     le nombre de points inferieurs total suivant z (ntz).
!     ntx =? nty =? ntz
!     hx represente le pas suivant x, hx = 1 / (ntx+1),
!     hy represente le pas suivant y, hy = 1 / (nty+1),
!     hz represente le pas suivant z, hz = 1 / (ntz+1).
!
!!!
!!!
!
!     Cree le          : Tue Sept 11 11:02:41 2001
!
!     Nature           : parallelisation par OpenMP
!
!     Auteur           : Guy Moebs (CRIHAN - France)
!                        <Guy.Moebs@crihan.fr>
!
!***********************************************************************!
!
      PROGRAM Poisson
!
      USE numerics
!
      IMPLICIT   NONE
!
INTERFACE
!
      INCLUDE 'i_initialisation.h'
!
      INCLUDE 'i_domaine.h'
!
      INCLUDE 'i_pmv.h'
!
      INCLUDE 'i_saxpy.h'
!
      INCLUDE 'i_prodscal.h'
!
      INCLUDE 'i_erreur.h'
!
END INTERFACE
!
!---Declaration  des  variables
!
!Solution  w
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  w
!
!Solution exacte
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  w_exact
!
!Second membre
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  b
!
!Vecteurs de travail
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  p
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  q
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  r
!
!Nombre  iterations  en  temps
      INTEGER   ::  it, itg
!
!Nombre  iterations  maximum  en  temps
      INTEGER   ::  it_max
!
!Nombre de threads
      INTEGER   ::  nbproc
!
!$    EXTERNAL      OMP_GET_NUM_THREADS
!$    INTEGER   ::  OMP_GET_NUM_THREADS
!
!Code erreur
      INTEGER   ::  ierr = 0
!
!Scalaires gradient conjugue
      REAL(rp)  ::  alpha_k
      REAL(rp)  ::  beta_k
      REAL(rp)  ::  rnorm_k
!
!Mesure  du  temps
      REAL(dp)  ::  t1, t2
      REAL(dp)  ::  timef
!
!Frequence des sorties
      INTEGER   ::  ifreq
!
!Consommation memoire
      REAL(rp)  ::  gmem
!
!Test  de  convergence
      LOGICAL   ::  convergence
      REAL(rp)  ::  prec
      REAL(rp)  ::  gerr
!
!Precision machine
      REAL(rp)  ::  machep
!
!***********************************************************************!
!
!Critere de convergence
      machep = EPSILON(machep)
      prec = MAX( machep, 1E-10_rp )
!
!Pour le processeur  0, connaitre  le  nombre  de  points  interieurs  total
      READ  (5,*) ntx, nty, ntz, rlambda, it_max, ifreq
!
! Affichage des parametres de l'execution
      WRITE (6,1010) 'Equation de Helmholtz 3D'
      WRITE (6,1020) 'Critere d arret, prec = ', prec
      WRITE (6,1020) 'Precision machine, machep = ', machep
      WRITE (6,1030) 'ntx    : nombre de points en x = ', ntx
      WRITE (6,1030) 'nty    : nombre de points en y = ', nty
      WRITE (6,1030) 'ntz    : nombre de points en z = ', ntz
      WRITE (6,1020) 'lambda : scalaire positif = ', rlambda
      WRITE (6,1030) 'it_max : nombre maximal d iterations pour la convergence = ', it_max
      WRITE (6,1030) 'ifreq  : frequence d affichage des residus = ', ifreq
!
 1010 FORMAT ( A, I5 )
 1020 FORMAT ( A, E12.4, / )
 1030 FORMAT (A, I6, /)
!
      IF ( rlambda < machep ) THEN
         WRITE (6,*) 'scalaire lambda positif strictement'
         STOP
      END IF
!
!Determiner  les  indices  du  domaine
      CALL domaine(  )
!
! Consommation memoire
      gmem = 0.046875_rp * REAL(ex-sx+3, rp) * REAL(ey-sy+3, rp) * REAL(ez-sz+3, rp)
      WRITE (6,1040) gmem
 1040 FORMAT ( 'Volume memoire de (w,b,p,q,r) : ', F14.0, ' ko', /)
!
!Allocation  dynamique  des  tableaux  u, u_exact, f
      ALLOCATE(  w      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 b      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 p      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 q      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 r      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 w_exact(sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), STAT = ierr )
      IF ( ierr /= 0 ) THEN
         WRITE (6,1010) 'ierr = ', ierr
         WRITE (6,1010) 'Probleme allocation'
         STOP
      END IF
!
!Mesure  du  temps  en  seconde  dans  la  boucle  en  temps
      t1  = timef ()
!
!Schema  iteratif  en  temps
      it  =  0
      gerr = zero
      rnorm_k = zero
      alpha_k = zero
!
!************************************************************************!
!
!Debut de la region parallele
!
!$OMP PARALLEL DEFAULT( NONE ) &
!$OMP SHARED( w, b, w_exact, p, q, r ) &
!$OMP SHARED( rnorm_k, alpha_k, beta_k, gerr ) &
!$OMP SHARED( itg, ifreq, nbproc, one, zero,  prec ) &
!$OMP PRIVATE( convergence ) &
!$OMP FIRSTPRIVATE( it, it_max )
!
!Initialisation  du  second  membre  et  de  la  solution  u  sur  les  bords
      CALL  initialisation(  w, p, r, b, w_exact, rnorm_k  )
!
!$OMP MASTER
         nbproc = 1
!$       nbproc = OMP_GET_NUM_THREADS ()
         WRITE (6,1050) 'Nombre de threads : ', nbproc
 1050 FORMAT (A, I5, /)
!
!Residus / erreurs
      WRITE (6,200)
      WRITE (6,201)
      WRITE (6,200)
      WRITE (6,202) 0, SQRT(rnorm_k)
  200 FORMAT ( '+-----------+------------+')
  201 FORMAT ( '| Iteration |   Residu   |')
  202 FORMAT ( '|  ', I6, '   |', E11.4, ' |')
!$OMP END MASTER
!
      convergence = ( SQRT(rnorm_k) < prec )
!
      DO  WHILE  ( (.NOT. convergence)  .AND.  (it  <  it_max) )
!
          it  =  it  +  1
!
!===================
! calcul de \alpha_k
!===================
!
! pmv   :  q_k = A.p_k
         CALL  pmv(  q,  p  )
!
! (q_k,p_k)
         CALL  prodscal( alpha_k,  q,  p  )
!
! alpha_k
!$OMP MASTER
         alpha_k = rnorm_k / alpha_k
         beta_k  = rnorm_k
         rnorm_k = zero
!$OMP END MASTER
!
!$OMP BARRIER
!
!================
! calcul de w_k+1
!================
         CALL  saxpy(  w,  one,  p,    alpha_k  )
!
!================
! calcul de r_k+1
!================
         CALL  saxpy(  r,  one,  q,  - alpha_k )
!
!==================
! calcul de \beta_k
!==================
!
! || r_k+1 ||^2
         CALL  prodscal(  rnorm_k,  r,  r  )
!
! \beta_k
!$OMP MASTER
         alpha_k = zero
         beta_k  = rnorm_k / beta_k
!$OMP END MASTER
!
!$OMP BARRIER
!
!================
! calcul de p_k+1
!================
         CALL  saxpy(  p,  beta_k,  r,  one  )
!
!Arret  du  programe  si  on  a  atteint  la  precision  souhaitee
         convergence  =  ( SQRT(rnorm_k)  <  prec )
!
!Affichage  du  residu  par  la  thread  maitresse
!$OMP MASTER
         itg = it
         IF  ( MOD(itg,ifreq)  ==  0 ) WRITE (6,202) itg, SQRT(rnorm_k)
!$OMP END MASTER
!
      END  DO
!
!************************************************************************!
!
!Validation avec une solution analytique
      CALL  erreur(  w,  w_exact,  gerr  )
!
!Fin   de la region parallele
!$OMP END PARALLEL
!
!Mesure  du  temps  a  la  sortie  de  la  boucle
      t2  =  ( timef () - t1 ) * 0.001_dp
!
!Affichage  des  performances
      IF ( MOD(itg, ifreq)  /=  0 )  THEN
         WRITE (6,202) itg, SQRT(rnorm_k)
      END IF
      WRITE (6,200)
!
      t1 = 15.0E-06_rp * REAL(ntx,rp) * REAL(nty,rp) * REAL(ntz,rp) * REAL(itg,rp) / t2
      WRITE (6,203) 'Erreur maximale : ', gerr
      WRITE (6,203) 'Temps ecoule : ', t2
      WRITE (6,204) 'Performances cumulees       : ', t1, ' MFlops'
      WRITE (6,204) 'Performances par processeur : ', t1 / REAL(nbproc,rp), ' MFlops'
!
  203 FORMAT ( A, E12.4 )
  204 FORMAT ( A, F8.2, A )
!
!Desallocation  des  tableaux
      DEALLOCATE(  w,  w_exact,  b,  q,  p,  r  )
!
      STOP
      END  PROGRAM  Poisson

