!>
!> @file ex6.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   Tranpsose of matrices partitionned on a 2d proc grid:
!     - A(n1/P1,n2,n3/P2)    -> AT(n1/P1,n3,n2/P2)
!     - B(n1,n2,n3/P1,n4/P2) -> BT(n1,n3,n2/P1,n4/P2)
!     - C(n1,n2,n3/P1,n4/P2) -> CT(n1,n4,n3/P1,n2/P2)
!     - D(n1,n2/P1,n3,n4/P2) -> DT(n1,n2/P1,n4,n3/P2)
!   n1, n2, n3, n4 NOT REQUIRED to be divided evenly by NPES
!
  USE pputils2
  USE futils
  IMPLICIT NONE
  INCLUDE "mpif.h"
  CHARACTER(len=32) :: file='ex4.h5'
  INTEGER :: fid
!
  INTEGER, PARAMETER :: ndims=2              ! N. of dims of proc. grid
  INTEGER :: ierr, me, npes
  INTEGER, DIMENSION(ndims) :: dims, coords
  LOGICAL :: periods(ndims), reorder
  INTEGER :: cart, cartcol, cartrow
!
  INTEGER ::  n1=8, n2=10, n3=6, n4=5, n1p, n2p, n3p, n4p, s1, s2, s3, s4
  DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: a3, a3t
  DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: b4, b4t, c4, c4t, d4, d4t
  DOUBLE PRECISION :: x
  INTEGER :: i, j, k, l, iglob, jglob, kglob, lglob, kerrors, nerrors
!================================================================================
!
!   Init MPI
  CALL mpi_init(ierr)
  CALL mpi_comm_size(MPI_COMM_WORLD, npes, ierr)
  CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
!
!   Create cartesian topololy
  dims    = (/4, 3/)
  periods = (/.FALSE., .TRUE./)
  reorder = .FALSE.
  IF( PRODUCT(dims) .NE. npes ) THEN
     IF( me .EQ. 0 ) THEN
        PRINT*,  PRODUCT(dims), " processors required!"
        CALL mpi_abort(MPI_COMM_WORLD, -1, ierr)
     END IF
  END IF
  CALL mpi_cart_create(MPI_COMM_WORLD, ndims, dims, periods, reorder, cart, ierr)
  CALL mpi_cart_coords(cart, me, ndims, coords, ierr)
  CALL mpi_cart_sub(cart, (/.TRUE., .FALSE. /), cartcol, ierr)
  CALL mpi_cart_sub(cart, (/.FALSE., .TRUE. /), cartrow, ierr)
!
!   Define local array A3
  CALL dist1d(cartcol, 0, n1, s1, n1p)
  CALL dist1d(cartrow, 0, n2, s2, n2p)
  CALL dist1d(cartrow, 0, n3, s3, n3p)
  ALLOCATE( a3(n1p,n2,n3p), a3t(n1p,n3,n2p) )
  a3 = 0
  a3t = 0
  DO i=1,n1p
     iglob = s1 + i
     DO j=1,n2
        DO k=1,n3p
           kglob = s3 + k
           a3(i,j,k) = 10000*iglob + 100*j + kglob
        END DO
     END DO
  END DO
  IF( me .EQ. 0 ) THEN
     WRITE(*,'(a,4i4)') 'Global dimension of matrix a', n1, n2, n3
  END IF
!
!   Tranpose A(n1/P1,n2,n3/P2) -> AT(n1/P1,n3,n2/P2)
  CALL pptransp(cartrow, a3, a3t, 2, 3)
!
  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
!
!   Check A3T
  kerrors = 0
  DO i=1,n1p
     iglob = s1 + i
     DO j=1,n2p
        jglob = s2 + j
        DO k=1,n3
           x = 10000*iglob + 100*jglob + k
           IF( x .NE. a3t(i,k,j) ) kerrors = kerrors+1
        END DO
     END DO
  END DO
  CALL mpi_reduce(kerrors, nerrors, 1, MPI_INTEGER, MPI_SUM, 0, &
       &          MPI_COMM_WORLD, ierr)
  IF( me .EQ. 0 ) WRITE(*,'(a,i6)') 'nerrors checking a3', nerrors
!
!   Define local array B4
  CALL dist1d(cartcol, 0, n2, s2, n2p)
  CALL dist1d(cartcol, 0, n3, s3, n3p)
  CALL dist1d(cartrow, 0, n4, s4, n4p)
  ALLOCATE( b4(n1,n2,n3p,n4p), b4t(n1,n3,n2p,n4p) )
  b4 = 0
  b4t = 0
  DO i=1,n1
     DO j=1,n2
        DO k=1,n3p
           kglob = s3 + k
           DO l=1,n4p
              lglob = s4 + l
              b4(i,j,k,l) = 1000000*i + 10000*j + 100*kglob +lglob
           END DO
        END DO
     END DO
  END DO
  IF( me .EQ. 0 ) THEN
     WRITE(*,'(a,4i4)') 'Global dimension of matrix b', n1, n2, n3, n4
  END IF
!
!   Tranpose B(n1,n2,n3/P1,n4/P2) -> BT(n1,n3,n2/P1,n4/P2)
  CALL pptransp(cartcol, b4, b4t, 2, 3)
!
  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
!
!   Check B4T
  kerrors = 0
  DO i=1,n1
     DO j=1,n2p
        jglob = s2 + j
        DO k=1,n3
           DO l=1,n4p
              lglob = s4 + l
              x = 1000000*i + 10000*jglob + 100*k + lglob
              IF( x .NE. b4t(i,k,j,l) ) kerrors = kerrors+1
           END DO
        END DO
     END DO
  END DO
  CALL mpi_reduce(kerrors, nerrors, 1, MPI_INTEGER, MPI_SUM, 0, &
       &          MPI_COMM_WORLD, ierr)
  IF( me .EQ. 0 ) WRITE(*,'(a,i6)') 'nerrors checking b4', nerrors
!
!   Define local array C4
  CALL dist1d(cartrow, 0, n2, s2, n2p)
  CALL dist1d(cartcol, 0, n3, s3, n3p)
  CALL dist1d(cartrow, 0, n4, s4, n4p)
  ALLOCATE( c4(n1,n2,n3p,n4p), c4t(n1,n4,n3p,n2p) )
  c4 = 0
  c4t = 0
  DO i=1,n1
     DO j=1,n2
        DO k=1,n3p
           kglob = s3 + k
           DO l=1,n4p
              lglob = s4 + l
              c4(i,j,k,l) = 1000000*i + 10000*j + 100*kglob +lglob
           END DO
        END DO
     END DO
  END DO
  IF( me.EQ. 0 ) THEN
     WRITE(*,'(a,4i4)') 'Global dimension of matrix c', n1, n2, n3, n4
  END IF
!
!   Tranpose C(n1,n2,n3/P1,n4/P2) -> CT(n1,n4,n3/P1,n2/P2)
  CALL pptransp(cartrow, c4, c4t, 2, 4)
!
  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
!
!   Check C4T
  kerrors = 0
  DO i=1,n1
     DO j=1,n2p
        jglob = s2 + j
        DO k=1,n3p
           kglob = s3 + k
           DO l=1,n4
              x = 1000000*i + 10000*jglob + 100*kglob + l
              IF( x .NE. c4t(i,l,k,j) ) kerrors = kerrors+1
           END DO
        END DO
     END DO
  END DO
  CALL mpi_reduce(kerrors, nerrors, 1, MPI_INTEGER, MPI_SUM, 0, &
       &          MPI_COMM_WORLD, ierr)
  IF( me .EQ. 0 ) WRITE(*,'(a,i6)') 'nerrors checking c4', nerrors
!
!   Define local array D4
  CALL dist1d(cartcol, 0, n2, s2, n2p)
  CALL dist1d(cartrow, 0, n3, s3, n3p)
  CALL dist1d(cartrow, 0, n4, s4, n4p)
  ALLOCATE( d4(n1,n2p,n3,n4p), d4t(n1,n2p,n4,n3p) )
  d4 = 0
  d4t = 0
  DO i=1,n1
     DO j=1,n2p
        jglob = s2 + j
        DO k=1,n3
           DO l=1,n4p
              lglob = s4 + l
              d4(i,j,k,l) = 1000000*i + 10000*jglob + 100*k +lglob
           END DO
        END DO
     END DO
  END DO
  IF( me.EQ. 0 ) THEN
     WRITE(*,'(a,4i4)') 'Global dimension of matrix d', n1, n2, n3, n4
  END IF
!
!   Tranpose D(n1,n2/P1,n3,n4/P2) -> DT(n1,n2/P1,n4,n3/P2)
  CALL pptransp(cartrow, d4, d4t, 3, 4)
!
  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
!
!   Check D4T
  kerrors = 0
  DO i=1,n1
     DO j=1,n2p
        jglob = s2 + j
        DO k=1,n3p
           kglob = s3 + k
           DO l=1,n4
              x = 1000000*i + 10000*jglob + 100*kglob + l
              IF( x .NE. d4t(i,j,l,k) ) kerrors = kerrors+1
           END DO
        END DO
     END DO
  END DO
  CALL mpi_reduce(kerrors, nerrors, 1, MPI_INTEGER, MPI_SUM, 0, &
       &          MPI_COMM_WORLD, ierr)
  IF( me .EQ. 0 ) WRITE(*,'(a,i6)') 'nerrors checking d4', nerrors
!
!   Write to file
!
  CALL creatf(file, fid, mpicomm=cart)
  CALL putarrnd(fid, '/a3' , a3, (/1,3/) )
  CALL putarrnd(fid, '/a3t', a3t,(/1,3/) )
  CALL putarrnd(fid, '/b4' , b4, (/3,4/) )
  CALL putarrnd(fid, '/b4t', b4t,(/3,4/) )
  CALL putarrnd(fid, '/c4' , c4, (/3,4/) )
  CALL putarrnd(fid, '/c4t', c4t,(/3,4/) )
  CALL putarrnd(fid, '/d4' , d4, (/2,4/) )
  CALL putarrnd(fid, '/d4t', d4t,(/2,4/) )
!
!   Clean up and quit
  DEALLOCATE(a3, a3t)
  DEALLOCATE(b4, b4t)
  DEALLOCATE(c4, c4t)
  DEALLOCATE(d4, d4t)
  CALL closef(fid)
  CALL mpi_finalize(ierr)
END PROGRAM main