!>
!> @file tp2p_mat.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
  USE pardiso_bsplines
  IMPLICIT NONE
  INCLUDE 'mpif.h'
!
  INTEGER :: me, npes, ierr
  INTEGER :: next
  INTEGER :: i, j, rank
  DOUBLE PRECISION :: val
  DOUBLE PRECISION, ALLOCATABLE :: arow(:)
  TYPE(pardiso_mat) :: mat
!
  CALL mpi_init(ierr)
  CALL mpi_comm_size(MPI_COMM_WORLD, npes, ierr)
  CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
  IF(npes.NE.2) THEN
     PRINT*, 'Should run with 2 procs!'
     CALL mpi_abort(MPI_COMM_WORLD, -1, ierr)
  END IF
!
!   Define local matrix
!
  rank = npes
  CALL init(npes, 0, mat)
  DO i=1,rank  ! Fill row me+1
     val = me+1
     j = me+1
     CALL updtmat(mat, i, me+1, val)
  END DO
!
!  Exchange matrix 
!
  CALL disp_mat('Original matrix')
  next=MODULO(me+1,2)
!
  IF(me.EQ.0) THEN
     CALL p2p_mat(mat, 1, 'send', 'put', MPI_COMM_WORLD)
  ELSE
     CALL p2p_mat(mat, 0, 'recv', 'put', MPI_COMM_WORLD)
  END IF
  CALL disp_mat('Matrix after 0->1/put')
!
  CALL p2p_mat(mat, next, 'sendrecv', 'put', MPI_COMM_WORLD)
  CALL disp_mat('Matrix after sendrev/put')
!
  CALL p2p_mat(mat, next, 'sendrecv', 'updt', MPI_COMM_WORLD)
  CALL disp_mat('Matrix after sendrev/updt')
!
  IF(me.EQ.1) THEN
     CALL p2p_mat(mat, 0, 'send', 'updt', MPI_COMM_WORLD)
  ELSE
     CALL p2p_mat(mat, 1, 'recv', 'updt', MPI_COMM_WORLD)
  END IF
  CALL disp_mat('Matrix after 1->0/updt')
!
  IF(me.EQ.1) THEN
     CALL p2p_mat(mat, 0, 'send', 'put', MPI_COMM_WORLD)
  ELSE
     CALL p2p_mat(mat, 1, 'recv', 'put', MPI_COMM_WORLD)
  END IF
  CALL disp_mat('Matrix after 1->0/put')
!
  CALL mpi_finalize(ierr)
CONTAINS
  SUBROUTINE disp_mat(str)
    CHARACTER(len=*), INTENT(in) :: str
    INTEGER :: p
    DO p=0,npes-1
       IF(me.EQ.p) THEN
          WRITE(*,'(a,i3.3)') str//' on PE', me
          CALL to_mat(mat, nlkeep=.TRUE.)
          ALLOCATE(arow(mat%rank))
          DO i=1,mat%rank
             CALL getrow(mat, i, arow)
             WRITE(*,'(10f8.2)') arow
          END DO
          DEALLOCATE(arow)
          CALL FLUSH(6)
       END IF
       CALL mpi_barrier(MPI_COMM_WORLD, ierr)
    END DO
  END SUBROUTINE disp_mat
END PROGRAM main