!>
!> @file poisson_mumps.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 mumps_bsplines
  USE cds
!
  IMPLICIT NONE
  INCLUDE 'mpif.h'
  TYPE(mumps_mat) :: amat
  TYPE(cds_mat) :: amat_cds
  DOUBLE PRECISION, ALLOCATABLE :: rhs(:), sol(:), arow(:)
  INTEGER :: nx=5, ny=4
  INTEGER :: n
  INTEGER :: i, j, irow
  INTEGER :: ierr, me
  INTEGER, ALLOCATABLE :: dists(:)
  DOUBLE PRECISION :: t0
!
  CALL mpi_init(ierr)
  CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
!
  IF(me.EQ.0) THEN
     WRITE(*,'(a)', advance='no') 'Enter nx, ny: '
     READ(*,*) nx, ny
  END IF
  CALL mpi_bcast(nx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
  CALL mpi_bcast(ny, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
  n = nx*ny  ! Rank of the matrix
  ALLOCATE(rhs(n))
  ALLOCATE(sol(n))
  ALLOCATE(arow(n))
!
  WRITE(*,'(/a)') 'Mumps using CSR mat ...'
  CALL init(n, 1, amat)
!
!  Construct the matrix and RHS
!
  t0 = mpi_wtime(0)
  DO j=1,ny
     DO i=1,nx
        arow = 0.0d0
        irow = numb(i,j)
        arow(irow) = 4.0d0
        IF(i.GT.1)  arow(numb(i-1,j)) = -1.0d0
        IF(i.LT.nx) arow(numb(i+1,j)) = -1.0d0
        IF(j.GT.1)  arow(numb(i,j-1)) = -1.0d0
        IF(j.LT.ny) arow(numb(i,j+1)) = -1.0d0
        CALL putrow(amat, irow, arow)
        rhs(irow) = SUM(arow)   ! => the exact solution is 1
     END DO
  END DO
!
  WRITE(*,'(a,i6)') 'Rank of matrix', n
  WRITE(*,'(a,i6)') 'Number of non-zeros of matrix', get_count(amat)
  WRITE(*,'(a,1pe12.3)') 'Matrix construction time (s)', mpi_wtime()-t0
!
!  Factor the amat matrix (Reordering, symbolic and numerical factorization)
!
  t0 = mpi_wtime(0)
  CALL factor(amat, nlmetis=.TRUE.)
  sol=rhs
  CALL bsolve(amat, sol)
  WRITE(*,'(a,1pe12.3)') 'Direct solve time (s)', mpi_wtime()-t0
!
  IF(me.EQ.0) THEN
     WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
  END IF
  CALL destroy(amat)
!
!   CDS matrix
!
  WRITE(*,'(/a)') 'Mumps using CDS mat ...'
  IF(ALLOCATED(dists)) DEALLOCATE(dists)
  ALLOCATE(dists(-2:2))
  dists = [-nx, -1, 0, 1, nx]
  WRITE(*,'(a/(20i4))') 'dists used in INIT =', dists  
  CALL init(n, dists, 1, amat_cds)
!
  t0 = mpi_wtime(0)
  DO j=1,ny
     DO i=1,nx
        arow = 0.0d0
        irow = numb(i,j)
        amat_cds%val(irow,0) = 4.0d0
        IF(i.GT.1)  amat_cds%val(irow,-1) = -1.0d0
        IF(i.LT.nx) amat_cds%val(irow,+1) = -1.0d0
        IF(j.GT.1)  amat_cds%val(irow,-2) = -1.0d0
        IF(j.LT.ny) amat_cds%val(irow,+2) = -1.0d0
     END DO
  END DO
  WRITE(*,'(a,1pe12.3)') 'Matrix construction time (s)', mpi_wtime()-t0
!
!   Compute dists of amat
  PRINT*, 'stat of mata%mat', ASSOCIATED(amat%mat)
  PRINT*, 'rank of mata', amat%mat%rank
  CALL mstruct(amat%mat, dists)
  WRITE(*,'(A/(20i4))') 'dists from MSTRUCT=', dists  
!
  t0 = mpi_wtime(0)
  CALL cds2mumps(amat_cds, amat)
  CALL factor(amat, debug=.FALSE.)
  sol = rhs
  CALL bsolve(amat, sol)
  WRITE(*,'(a,1pe12.3)') 'Direct solve time (s)', mpi_wtime()-t0
!
  IF(me.EQ.0) THEN
     WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
  END IF
!
!   Clean up
!
  DEALLOCATE(rhs)
  DEALLOCATE(sol)
  DEALLOCATE(arow)
  CALL destroy(amat)
  CALL mpi_finalize(ierr)
CONTAINS
  SUBROUTINE mstruct(mat, dists)
    TYPE(spmat), INTENT(in) :: mat
    INTEGER, ALLOCATABLE, INTENT(inout) :: dists(:)
    TYPE(elt), POINTER :: t
    INTEGER :: n, i, j0
    j0 = LBOUND(dists,1)
    n = mat%rank
    PRINT*, 'rank of mat', n
    DO i=1,n                    ! scan the matrix rows   
       t => mat%row(i)%row0
       DO WHILE(ASSOCIATED(t))  ! walk thru the linked list row(i)
          j = t%index
          IF(ABS(t%val) .LE. EPSILON(0.0d0)) THEN
             dists(j0) = t%index-i ! distance from main diag
             j0 = j0+1
          END IF
          t => t%next
       END DO
    END DO
  END SUBROUTINE mstruct
  INTEGER FUNCTION numb(i,j)
!
!  One-dimensional numbering
!  Number first x then y
!
    INTEGER, INTENT(in) :: i, j
    numb = (j-1)*nx + i
  END FUNCTION numb
END PROGRAM main