!************************************************************************!
!Resolution de l'equation de la chaleur sur le domaine [0,1]x[0,1]x[0,1]
!     discretisation spatiale : differences finies 2nd ordre
!     discretisation temporelle : Crank-Nicholson, 2nd ordre
!     avec comme solveur le Gradient Conjugue.
!     d(u)/dt - \lambda * Delta  u =  f(x,y,z)
!     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       : Mon Sept 17  9:32:41 2001
!
!     Nature           : parallelisation a deux niveaux
!                        decomposition de domaines avec MPI
!                        partage du travail        avec OpenMP
!                        version 2 : decoupage des boucles a la main
!                        passage du code en 3D
!                        conversion Poisson -> Chaleur
!                        conversion jacobi  -> gradient conjugue
!
!     Auteur           : Guy Moebs (CRIHAN - France)
!                        <Guy.Moebs@crihan.fr>
!
!***********************************************************************!
!
      PROGRAM  Chaleur
!
      USE MPI
!
      USE numerics
!
      IMPLICIT   NONE
!
      INCLUDE 'comthd.h'
!
INTERFACE
!
      INCLUDE     'i_communication.h'
!
      INCLUDE     'i_cree_type.h'
!
      INCLUDE     'i_domaine.h'
!
      INCLUDE     'i_fctfx.h'
!
      INCLUDE     'i_scdmb.h'
!
      INCLUDE     'i_gradconj.h'
!
      INCLUDE     'i_nrmerr.h'
!
      INCLUDE     'i_voisinage.h'
!
END INTERFACE
!
!---Declaration  des  variables
!
!Solution  u
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  u
!
!Solution  exacte
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  u_exact
!
!Second membre
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  b
!
!Forcage  exterieur
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  f1
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  f2
!
!Vecteurs de travail
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  p
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  q
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  r
!
! Tableaux des repartitions
      INTEGER,  DIMENSION(:),  ALLOCATABLE  ::  npts
      INTEGER,  DIMENSION(:),  ALLOCATABLE  ::  ntproc
!
!Nombre  iterations  en  temps
      INTEGER   ::  itg
!
!Nombre  d'iterations  maximum  gradient  conjugue
      INTEGER   ::  it_max
!
!Nombre  d'iterations  en  temps
      INTEGER   ::  nbiter
!
!Nombre  de  threads
      INTEGER   ::  ir, iq
      INTEGER   ::  iproc
!
!$    EXTERNAL      OMP_GET_NUM_THREADS
!$    INTEGER   ::  OMP_GET_NUM_THREADS
!
!$    EXTERNAL      OMP_GET_THREAD_NUM
!$    INTEGER   ::  OMP_GET_THREAD_NUM
!
!Scalaires  gradient  conjugue
      REAL(rp)  ::  rnorm_k
      INTEGER   ::  itgc
      INTEGER   ::  nbitgc
!
!Scalaire  Crank  Nicholson
      REAL(rp)  ::  rnudt2
!
!Normes  diverses
      REAL(rp),  DIMENSION(2)  ::  gnrm2
      REAL(rp),  DIMENSION(2)  ::  gnrmoo
!
!Temps  courant
      REAL(rp)  ::  tps
!
!Frequence  des  sorties
      INTEGER   ::  ifreq
!
!Consommation  memoire
      REAL(rp)  ::  gmem
!
!Convergence  gradient  conjugue
      REAL(rp)  ::  prec
!
!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
!
!Compteur
      INTEGER   ::  j
!
!Mesure  du  temps
      REAL(rp)  ::  t1, t2
!
!Volume memoire
      REAL(rp)  ::  ilmem
!
!Test  de  convergence
      LOGICAL   ::  convergence
!
!Constantes  MPI
!
      INTEGER   ::  imyrank
      INTEGER   ::  comm3d
      INTEGER   ::  ip_xy,  ip_xz,  ip_yz
      INTEGER,  DIMENSION(3)  ::  dims
      LOGICAL,  DIMENSION(3)  ::  periods
!
     !***********************************************************************!
!
!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  )
!
!Lecture des parametres
      machep = MAX( EPSILON(machep), 1E-15_rp )
!
      IF  ( myrank  ==  0 )  THEN
         OPEN(19,  FILE='C3D.dat',  STATUS='OLD',  FORM='FORMATTED')
         READ(19, 101) nbthd
         READ(19, 101) ntx
         READ(19, 101) nty
         READ(19, 101) ntz
         READ(19, 101) nbiter
         READ(19, 102) dt
         READ(19, 102) rlambda
         READ(19, 102) prec
         READ(19, 101) it_max
         READ(19, 101) ifreq
         CLOSE(19)
!
  101 FORMAT ( I5 )
  102 FORMAT ( E14.4 )
!
         WRITE (6,201) 'Equation de la chaleur 3D'
         WRITE (6,202) 'machep = ', machep
         WRITE (6,201) 'Nombre de threads OpenMP : ', nbthd
         WRITE (6,201) 'Nombre de points en X : ', ntx
         WRITE (6,201) 'Nombre de points en Y : ', nty
         WRITE (6,201) 'Nombre de points en Z : ', ntz
!
         WRITE (6,201) 'Nombre d iterations en temps : ', nbiter
         WRITE (6,202) 'Pas de temps : ', dt
         WRITE (6,202) 'Viscosite cinematique : ', rlambda
!
         WRITE (6,202) 'Precision gradient    : ', prec
         WRITE (6,201) 'Nbre max d iterations : ', it_max
         WRITE (6,201) 'Frequence des sorties : ', ifreq
!
  201 FORMAT ( A, I5 )
  202 FORMAT ( A, E12.4 )
!
         WRITE (6,1010) 'Nombre de processus MPI : ', nbproc
!
         IF ( rlambda < machep ) THEN
              WRITE (6,*) 'scalaire lambda positif strictement'
              CALL MPI_ABORT ( MPI_COMM_WORLD, 1, ierr )
              STOP
          END IF
!
      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(  nbthd,    1,  MPI_INTEGER,   0,  MPI_COMM_WORLD,  ierr  )
      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(  nbiter,   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  )
      CALL  MPI_BCAST(  dt,       1,  ITYPE_REAL,    0,  MPI_COMM_WORLD,  ierr  )
      CALL  MPI_BCAST(  rlambda,  1,  ITYPE_REAL,    0,  MPI_COMM_WORLD,  ierr  )
      CALL  MPI_BCAST(  prec,     1,  ITYPE_REAL,    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  )
      write (6,*) 'dims = ', dims
!
!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.0625_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), &
                 b  (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 f1 (sx  :ex  , sy  :ey  , sz  :ez  ), &
                 f2 (sx  :ex  , sy  :ey  , sz  :ez  ), &
                 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), &
                 u_exact(sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 npts   (0:nbthd-1), ntproc (0:nbthd),     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
!
!Repartition des iterations selon z
      iq = (ez - sz + 1) / nbthd
!
      DO iproc = 0, nbthd - 1
         npts(iproc) = iq
      END DO
!
      ir = (ez - sz + 1) - nbthd * iq
      IF ( ir /= 0 ) THEN
         DO iproc = 0, ir - 1
            npts(iproc) = npts(iproc) + 1
         END DO
      END IF
!
      ntproc (0) = sz
      DO iproc = 0, nbthd - 1
         ntproc(iproc+1) = ntproc(iproc) + npts(iproc)
      END DO
!
      rnudt2 = 0.5_rp * rlambda * dt
!
      itgc   = 0
      nbitgc = 0
!
      rnorm_k = zero
      tps     = zero
!
!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  )
!
!************************************************************************!
!
!Debut de la region parallele
!
      CALL OMP_SET_NUM_THREADS ( nbthd )
!
!$OMP PARALLEL DEFAULT( NONE ) &
!$OMP SHARED( u, b, u_exact, f1, f2, p, q, r ) &
!$OMP SHARED( tps, dt, rnorm_k, itgc, nbitgc ) &
!$OMP SHARED( gnrm2,  gnrmoo ) &
!$OMP SHARED( nbproc, ntproc, nbthd, nbiter, prec, machep ) &
!$OMP SHARED( ifreq, it_max, rnudt2, rlambda ) &
!$OMP SHARED( ip_xy, ip_xz, ip_yz, voisin, comm3d ) &
!$OMP SHARED( t2, t1, myrank, ierr ) &
!$OMP PRIVATE( itg )
!
!$OMP MASTER
         nbthd = 1
!$       nbthd = OMP_GET_NUM_THREADS ()
         WRITE (6,1040) 'Processus MPI : ', myrank, ' Nombre de threads OpenMP : ', nbthd
!$OMP END MASTER
 1040 FORMAT ( 2(A, I5), /)
!
      mythd = 0
!$    mythd = OMP_GET_THREAD_NUM ()
!
      kmin = ntproc(mythd)
      kmax = ntproc(mythd+1) - 1
!
!Initialisation  des  matrices
      CALL  fctfx(  u,  f1,  f2,  p,  u_exact,  tps  )
!
!Erreurs initiales
      CALL  nrmerr(  gnrm2,  gnrmoo,  u,  u_exact,  tps,  comm3d  )
!
      IF ( myrank == 0 ) THEN
!$OMP MASTER
         WRITE (6, 401) 'TPS', 'uL2', 'uLoo', 'eL2', 'eLoo', 'itGC', 'res'
         WRITE (6, 402) tps,  gnrm2(1),  gnrmoo(1),  gnrm2(2),  gnrmoo(2),  itgc,  SQRT(rnorm_k)
!$OMP END MASTER
  401 FORMAT ( 4x, A3, 7x, A3, 8x, A4, 9x, A3, 8x, A4, 7x, A4, 5x, A3 )
  402 FORMAT ( F8.4, 4E12.4, 4x, I4, E12.4 )

      END IF
!
!Mesure  du  temps  en  seconde  dans  la  boucle  en  temps
!$OMP MASTER
      t1  =  MPI_WTIME()
!$OMP END MASTER
!
!=================================
! Debut  de  la  boucle  en  temps
!=================================
!
      DO itg = 1, nbiter
!
!Echange  des  points  aux  interfaces  pour  u
!$OMP MASTER
         CALL communication(  u,  ip_xy,  ip_xz,  ip_yz,  voisin,  comm3d  )
!
         nbitgc = nbitgc + itgc
         itgc = 0
         rnorm_k = zero
!
!$OMP END MASTER
!
!$OMP BARRIER
!
!Calcul  second  membre  Crank  Nicholson
         CALL  scdmb(  b,  f1,  f2,  u,  tps,  rnudt2  )
!
!Avancement en temps
!$OMP MASTER
         tps = tps + dt
!$OMP END MASTER
!
!Resolution  gradient  conjugue
         CALL  gradconj(  u,  b,  p,  q,  r,  it_max,  prec,  itgc,  rnorm_k, &
                          rnudt2,  ip_xy,  ip_xz,  ip_yz,  voisin,  comm3d  )
!
!Calcul des normes
         IF  ( MOD(itg,ifreq)  ==  0 ) THEN
!
            CALL  nrmerr(  gnrm2,  gnrmoo,  u,  u_exact,  tps,  comm3d  )
!
            IF ( myrank  ==  0 ) THEN
!$OMP MASTER
               WRITE (6, 402) tps,  gnrm2(1),  gnrmoo(1),  gnrm2(2),  gnrmoo(2),  itgc,  SQRT(rnorm_k)
!$OMP END MASTER
!
            END IF
!
         END IF
!
!=================================
! Fin    de  la  boucle  en  temps
!=================================
!
      END  DO
!
!Mesure  du  temps  a  la  sortie  de  la  boucle
!$OMP MASTER
      t2  =  MPI_WTIME() - t1
!$OMP END MASTER
!
!************************************************************************!
!
!Fin   de la region parallele
!$OMP END PARALLEL
!
         nbitgc = nbitgc + itgc
!
!Affichage  du  temps  de  convergence  par  le  processeur  0
      IF  ( myrank  ==  0 )  THEN
         t1 = REAL(ntx,rp) * REAL(nty,rp) * REAL(ntz,rp) * ( &
              23.0E-06_rp  * REAL(nbitgc,rp) &                  ! gradconj
            + 30.0E-06_rp  * REAL(nbiter)    &                  ! scdmb
            + 12.0E-06_rp  * REAL(nbiter/ifreq + 1, rp) ) / t2  ! nrmerr
         WRITE (6,202) 'Temps ecoule : ', t2
         WRITE (6,203) 'Performances cumulees             : ', t1, ' MFlops'
         WRITE (6,203) 'Performances par processus-thread : ', t1 / REAL(nbproc*nbthd,rp), ' MFlops'
!
  203 FORMAT ( A, F8.2, A )
!
      END  IF
!
!Desallocation  des  tableaux
      DEALLOCATE(  u,  u_exact,  b,  f1,  f2,  p,  q,  r,  igmem,  npts,  ntproc  )
!
!Liberation  du  type  type_ligne  et  du  communicateur  comm2d
       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)
!
      END  PROGRAM  Chaleur
