!************************************************************************!
!Resolution de l'equation de Poisson en utilisant Jacobi
!     sur le domaine [0,1]x[0,1]x[0,1] par une methode aux differences finies,
!     avec comme solveur Jacobi.
!     - 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).
!
!     Pour chaque processeur :
!     1) decomposer le domaine
!     2) connaitre ses 6 voisins
!     3) echanger les points aux interfaces
!     4) calculer
!
!     Auteur           : Isabelle DUPAYS (CNRS/IDRIS - France)
!                        <Isabelle.Dupays@idris.fr>
!     Cree le          : Fri nov 15 10:32:45 1996
!
!!!
!!!
!
!     Modifie le       : Tue May 15 18:02:41 2001
!
!     Nature           : passage du code en 3D
!
!     Auteur           : Guy Moebs (CRIHAN - France)
!                        <Guy.Moebs@crihan.fr>
!
!***********************************************************************!
!
      PROGRAM Poisson
!
      USE MPI
!
      USE numerics
!
      IMPLICIT   NONE
!
INTERFACE
!
      INCLUDE 'i_voisinage.h'
!
      INCLUDE 'i_domaine.h'
!
      INCLUDE 'i_cree_type.h'
!
      INCLUDE 'i_initialisation.h'
!
      INCLUDE 'i_communication.h'
!
      INCLUDE 'i_calcul.h'
!
      INCLUDE 'i_converge.h'
!
      INCLUDE 'i_erreur.h'
!
END INTERFACE
!
!---Declaration  des  variables
!
!Solution  u  a  l'iteration  n  et  n+1
      REAL(rp),  DIMENSION(:,:,:,:),  ALLOCATABLE  ::  u
!
!Solution  exacte
      REAL(rp),  DIMENSION(:,:,:),    ALLOCATABLE  ::  u_exact
!
!Second membre
      REAL(rp),  DIMENSION(:,:,:),    ALLOCATABLE  ::  f
!
!Consommation memoire
      REAL(rp),  DIMENSION(:),        ALLOCATABLE  ::  igmem
!
!Tableau  contenant  les  voisins  du  domaine  courant
      INTEGER,  PARAMETER                  ::  NB_VOISINS  =  6
!
      INTEGER,  DIMENSION(NB_VOISINS)      ::  voisin
      INTEGER,  PARAMETER                  ::  N=1,  E=2,   S=3
      INTEGER,  PARAMETER                  ::  W=4,  AV=5,  AR=6
!
!Nombre  iterations  en  temps
      INTEGER   ::  it
!
!Nombre  iterations  maximum  en  temps
      INTEGER   ::  it_max
!
!Frequence des sorties
      INTEGER   ::  ifreq
!
!Convergence  et  erreur
      REAL(rp)  ::  gdiff
      REAL(rp)  ::  gerr
!
!Mesure  du  temps
      REAL(rp)  ::  t1, t2
!
!Compteur
      INTEGER   ::  j
!
!Volume memoire
      REAL(rp)  ::  ilmem
!
!Test  de  convergence
      LOGICAL   ::  convergence
      REAL(rp)  ::  prec
!
!Constantes  MPI
!
      INTEGER   ::  imyrank
      INTEGER   ::  comm3d
      INTEGER   ::  ip_xy,  ip_xz,  ip_yz
      INTEGER,  DIMENSION(3)  ::  dims
      LOGICAL,  DIMENSION(3)  ::  periods
!
!Precision machine
      REAL(rp)  ::  machep
!
     !***********************************************************************!
!
!Initialisation  de  MPI
      CALL  MPI_INIT(ierr)
!
!Savoir  quel  processeur  je  suis
      CALL  MPI_COMM_RANK(  MPI_COMM_WORLD,  myrank,  ierr)
!
!Connaitre le nombre total de processeurs
      CALL  MPI_COMM_SIZE(  MPI_COMM_WORLD,  nbproc,  ierr)
!
!Critere de convergence
      machep = EPSILON(machep)
      prec = MAX( machep, 1E-10_rp )
!
!Pour le processeur  0, connaitre  le  nombre  de  points  interieurs  total
      IF (myrank  ==  0)  THEN
         READ  (5,*) ntx, nty, ntz, it_max, ifreq
         WRITE (6,1020) 'Equation de Poisson 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,1030) 'it_max : nombre maximal d iterations pour la convergence = ', it_max
         WRITE (6,1030) 'ifreq  : frequence d affichage des residus = ', ifreq
!
 1020 FORMAT ( A, E12.4, / )
 1030 FORMAT (A, I6, /)
      END IF
!
!Type de donnees
      IF ( rp == dp ) THEN
         ITYPE_REAL = MPI_DOUBLE_PRECISION
      ELSE
         ITYPE_REAL = MPI_REAL
      END IF
!
!Envoi de la valeur de ntx a tous les autres processeurs
      CALL  MPI_BCAST(  ntx,    1,  MPI_INTEGER,   0,  MPI_COMM_WORLD,  ierr  )
      CALL  MPI_BCAST(  nty,    1,  MPI_INTEGER,   0,  MPI_COMM_WORLD,  ierr  )
      CALL  MPI_BCAST(  ntz,    1,  MPI_INTEGER,   0,  MPI_COMM_WORLD,  ierr  )
      CALL  MPI_BCAST(  it_max, 1,  MPI_INTEGER,   0,  MPI_COMM_WORLD,  ierr  )
      CALL  MPI_BCAST(  ifreq,  1,  MPI_INTEGER,   0,  MPI_COMM_WORLD,  ierr  )
!
!Connaitre
!  le  nombre  de  processeurs  selon  x
!  le  nombre  de  processeurs  selon  y
!  le  nombre  de  processeurs  selon  z  en  fonction  du  nombre  total  de  processeurs
      dims(:)  =  0
      IF ( ntx == 1 ) dims(1) = 1
      IF ( nty == 1 ) dims(2) = 1
      IF ( ntz == 1 ) dims(3) = 1
      CALL  MPI_DIMS_CREATE(  nbproc,  3,  dims,  ierr)
!
!Creation  de  la  grille  de  processeurs  3D  sans  periodicite
      periods(:)  =  .FALSE.
!
      CALL  MPI_CART_CREATE( MPI_COMM_WORLD,  3,  dims,  periods,  .TRUE., &
                             comm3d,  ierr)
!
!Recherche  de  ses  6  voisins  pour  chaque  processeur
      CALL  voisinage(  comm3d,  voisin  )
!
!Determiner  les  indices  de  chaque  sous-domaine
      CALL domaine(  comm3d  )
!
! Consommation memoire en ko
      ilmem = 0.0234375_rp * REAL(ex-sx+3, rp) * REAL(ey-sy+3, rp) * REAL(ez-sz+3, rp)
      ALLOCATE(  igmem (0:nbproc-1), STAT=ierr)
      igmem (0:nbproc-1) = 0
!
      CALL  MPI_GATHER( ilmem, 1, ITYPE_REAL, &
                        igmem, 1, ITYPE_REAL, &
                        0, comm3d, ierr )
!
     CALL  MPI_COMM_RANK(  comm3d,  imyrank,  ierr)
      IF ( imyrank == 0 ) THEN
         WRITE (6,300)
         WRITE (6,301)
         WRITE (6,300)
         DO j = 0, nbproc - 1
            WRITE (6,302) j, igmem(j)
         END DO
         WRITE (6,300)
         WRITE (6,303) SUM( igmem(0:nbproc-1) )
         WRITE (6,300)
         WRITE (6,*)
  300 FORMAT ( '+-----------+------------+')
  301 FORMAT ( '| Processus | Mem. (ko)  |')
  302 FORMAT ( '|  ', I6, '   |', F12.1, '|')
  303 FORMAT ( '|   TOTAL   |', F12.1, '|')
      END IF
!
!Allocation  dynamique  des  tableaux
      ALLOCATE(  u      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1,  0:1), &
                 f      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 u_exact(sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), STAT = ierr )
      IF ( ierr /= 0 ) THEN
         WRITE (6,1010) 'rang = ', myrank
         WRITE (6,1010) 'ierr = ', ierr
         WRITE (6,1010) 'Probleme allocation'
         CALL MPI_ABORT ( MPI_COMM_WORLD, 1, ierr )
 1010 FORMAT ( A, I6 )
         STOP
      END IF
!
!Creation  du  type  type_ligne  pour  echanger  les  points  a  l'ouest  et  a  l'est
      CALL cree_type(  ip_xy,  ip_xz,  ip_yz  )
!
!Initialisation  du  second  membre  et  de  la  solution  u  sur  les  bords
      CALL  initialisation(  u,  f,  u_exact  )
!
      IF ( myrank == 0 ) THEN
         WRITE (6,200)
         WRITE (6,201)
         WRITE (6,200)
  200 FORMAT ( '+-----------+------------+')
  201 FORMAT ( '| Iteration |  EcartLoo  |')
      END IF
!
!Barriere
      CALL  MPI_BARRIER(  comm3d,  ierr  )
!
!Schema  iteratif  en  temps
      it  =  0
      convergence  =  .FALSE.
      gdiff = zero
!
!Mesure  du  temps  en  seconde  dans  la  boucle  iterative
      t1  =  MPI_WTIME()
!
      DO  WHILE  ( ( .NOT. convergence )  .AND.  (it  <  it_max)  .AND.  &
                   ( gdiff  <  50.0_rp ) )
!
         it  =  it  +  1
         j = iold
         iold = inew
         inew = j
!
!Echange  des  points  aux  interfaces  pour  u  a  l'iteration  it
         CALL communication(  u,  ip_xy,  ip_xz,  ip_yz,  voisin,  comm3d)
!
!Calcul  de  u  a  l'iteration  it+1
         CALL  calcul(  u,  f  )
!
!Calcul  de  l'ecart  sur  le  sous-domaine
         CALL  converge(  u,  gdiff,  comm3d  )
!
!Arret  du  programe  si  on  a  atteint  la  precision  souhaitee
         convergence  =  ( gdiff  <  prec )
!
!Affichage  pour  le  processeur  0  de  la  difference
         IF  ( (myrank  ==  0) .AND. (MOD(it, ifreq)  ==  0) )  THEN
            WRITE (6,202) it, gdiff
  202 FORMAT ( '|  ', I6, '   |', E11.4, ' |')
         END  IF
!
      END  DO
!
!Validation avec une solution analytique
      CALL  erreur( u,  u_exact,  gerr,  comm3d  )
!
!Mesure  du  temps  a  la  sortie  de  la  boucle
      t2  =  MPI_WTIME() - t1
!
!Affichage  du  temps  de  convergence  par  le  processeur  0
      IF  ( myrank  ==  0 )  THEN
         IF ( MOD(it, ifreq)  /=  0 )  THEN
            WRITE (6,202) it, gdiff
         END IF
         WRITE (6,200)
!
         t1 = 12.0E-06_rp * REAL(ntx,rp) * REAL(nty,rp) * REAL(ntz,rp) * REAL(it,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 )
      END  IF
!
!Desallocation  des  tableaux
      DEALLOCATE(  u,  f,  u_exact,  igmem  )
!
!Liberation  du  type  type_ligne  et  du  communicateur  comm3d
       CALL  MPI_TYPE_FREE( ip_xy,   ierr )
       CALL  MPI_TYPE_FREE( ip_xz,   ierr )
       CALL  MPI_TYPE_FREE( ip_yz,   ierr )
       CALL  MPI_COMM_FREE( comm3d,  ierr )
!
!Desactivation  de  MPI
      CALL  MPI_FINALIZE(ierr)
!
      STOP
      END  PROGRAM  Poisson

