!>
!> @file pde1d_eig_zmumps.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
!
!  Solving the following 1d differential  eqation using splines:
!
!  Solve the standard eigenvalue:
!      A*x = \lambda *x or inv(A)*x = 1/\lambda * x using Arpack and MUMPS.
!   where A is obtained from discretozation of
!     -d/dr[r d/dr] f = k^2 r^(k-1)
!
  USE bsplines
  USE mumps_bsplines
  USE futils
  IMPLICIT NONE
  INTEGER :: nx, nidbas, ngauss, kdiff
  INTEGER :: i, nrank
  LOGICAL :: nlppform
  DOUBLE PRECISION :: coefs(5)
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid
  TYPE(spline1d) :: splx
  TYPE(zmumps_mat) :: mat
!
  INTEGER :: ierr
  INTEGER :: fid
  CHARACTER(len=32) :: str
  DOUBLE PRECISION :: seconds, t0, tmat, tfac, tarpack
!
!   Arpack:  Solve the standard eigenvalue problem
!
  INTEGER :: nev = 10, ncv = 10
  LOGICAL :: nlinv = .FALSE.  ! Solve inv(A) = 1/\lambda * x if nlinv=.TRUE.
  CHARACTER(len=2) :: which='SM'
  INTEGER :: info=0   ! Use random vector to start the Arnoldi iterations
  INTEGER :: ido=0    ! Reverse communications
  LOGICAL :: rvec
  LOGICAL, ALLOCATABLE :: select(:)
  INTEGER :: iparam(11), ipntr(14), nconv
  DOUBLE PRECISION :: tol=0.0d0
  CHARACTER(len=1) :: bmat='I'
!
  INTEGER :: lworkl
  DOUBLE COMPLEX, ALLOCATABLE :: workl(:), workd(:), workev(:)
  DOUBLE COMPLEX, ALLOCATABLE :: eig_vals(:), eig_vecs(:,:), resid(:)
  DOUBLE COMPLEX :: sigma
  DOUBLE PRECISION, ALLOCATABLE :: rwork(:)
!
  INTERFACE
     SUBROUTINE dismat(spl, mat)
       USE bsplines
       USE mumps_bsplines
       TYPE(spline1d), INTENT(in) :: spl
       TYPE(zmumps_mat), INTENT(inout) :: mat
     END SUBROUTINE dismat
     SUBROUTINE disrhs(kdiff, spl, rhs)
       USE bsplines
       INTEGER, INTENT(in) :: kdiff
       TYPE(spline1d), INTENT(in) :: spl
       DOUBLE PRECISION, INTENT(out) :: rhs(:)
     END SUBROUTINE disrhs
     SUBROUTINE meshdist(coefs, x, nx)
       DOUBLE PRECISION, INTENT(in) :: coefs(5)
       INTEGER, INTENT(iN) :: nx
       DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
     END SUBROUTINE meshdist
  END INTERFACE
!
  NAMELIST /newrun/ nx, nidbas, ngauss, kdiff, nlppform, coefs, &
       &            nev, ncv, nlinv, which, tol
!===========================================================================
!              1.0 Prologue
!
  CALL mpi_init(ierr)
!
!   Read in data specific to run
!
  nx = 8     ! Number oh intevals in x
  nidbas = 3 ! Degree of splines
  ngauss = 4 ! Number of Gauss points/interval
  kdiff = 2  ! Exponent of differential problem
  nlppform = .TRUE. ! Use PPFORM for gridval or not
  coefs(1:5) = (/0.2d0, 1.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
                                                 ! see function FDIST in MESHDIST
!
  READ(*,newrun)
  WRITE(*,newrun)
!
!   Define grid on x axis
!
  ALLOCATE(xgrid(0:nx))
  xgrid(0) = 0.0d0
  xgrid(nx) = 1.0d0
  CALL meshdist(coefs, xgrid, nx)
  WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
  WRITE(*,'(a/(10f8.3))') 'Diff XGRID', (xgrid(i)-xgrid(i-1), i=1,nx)
  WRITE(*,'(a/(10f8.3))') 'Diff XGRID**2', (xgrid(i)**2-xgrid(i-1)**2, i=1,nx)
!===========================================================================
!              2.0 Discretize the PDE
!
!   Set up spline
!
  t0 = seconds()
  CALL set_spline(nidbas, ngauss, xgrid, splx, nlppform=nlppform)
  CALL get_dim(splx, nrank)       ! Rank of the FE matrix
  WRITE(*,'(a,l6)') 'splx%nlequid', splx%nlequid
!
!   FE matrix assembly
!
  WRITE(*,'(a,i6)') ' nrank', nrank
  CALL init(nrank, 1, mat)
  CALL dismat(splx, mat)
  CALL to_mat(mat)
  tmat = seconds() - t0
!
  mat%mumps_par%IRN => mat%mumps_par%IRN_loc ! Work around for single proc.
  PRINT*, 'nnz_loc', mat%nnz_loc
  PRINT*, 'mat%mumps_par%N', mat%mumps_par%N
  PRINT*, 'mat%mumps_par%NZ_loc', mat%mumps_par%NZ_loc
  PRINT*, 'size of mat%mumps_par%IRN', SIZE(mat%mumps_par%IRN)
  PRINT*, 'mat%istart,mat%iend', mat%istart,mat%iend
!  
  CALL creatf('pde1d_eig_zmumps.h5', fid, 'PDE1D Result File', real_prec='d')
  PRINT*, 'rank', mat%rank
  PRINT*, 'nnz', mat%nnz
!
  CALL putmat(fid,'/MAT',mat,'FE matrix')
!
  IF(nlinv) THEN
     t0 = seconds()
     CALL factor(mat)
     tfac = seconds()-t0
  END IF
!===========================================================================
!              3.0 Solve the standard eigenvalue problem
!
  lworkl = 3*ncv**2 + 5*ncv
  ALLOCATE(workl(lworkl))
  ALLOCATE(workd(3*nrank))
  ALLOCATE(workev(2*ncv))
  ALLOCATE(eig_vals(ncv), eig_vecs(nrank,ncv))
  ALLOCATE(resid(nrank))
  ALLOCATE(rwork(ncv))
!
  iparam(1) = 1   ! shfts
  iparam(3) = 300 ! Max. number of iterations
  iparam(7) = 1   ! Regular mode
!
! The reverse communication loop
!
  t0 = seconds()
  DO
     CALL znaupd  (ido, bmat, nrank, which, nev, tol, resid, ncv, &
          &        eig_vecs, nrank, iparam, ipntr, workd, workl, lworkl, &
          &        rwork, info )
!
     IF(ido .EQ. -1 .OR. ido .EQ. 1) THEN   ! Compute A*v
        CALL av(nrank, workd(ipntr(1)), workd(ipntr(2)))
        CYCLE
     END IF
!
     IF(info .LT. 0) THEN  ! Error
        PRINT*, 'Error in _naupd with info =', info
        EXIT
     ELSE
        rvec = .TRUE.
        ALLOCATE(select(ncv))
        CALL zneupd  (rvec, 'A', select, eig_vals, eig_vecs, nrank, &
             &       sigma, workev, bmat, nrank, which, nev, tol, resid, &
             &       ncv, eig_vecs, nrank, iparam, ipntr, workd, workl, lworkl,&
             &       rwork, ierr)
        IF(ierr .NE. 0) THEN
           PRINT*, 'Error in _neupd with ierr =', ierr
           EXIT
        ELSE
           nconv = iparam(5)
           PRINT*,'Number of converged eigenvalues', nconv
           IF(nlinv) THEN
              eig_vals(1:nconv) = (1.d0,0.0d0) / eig_vals(1:nconv)
           END IF
           WRITE(*,'(2(1pe12.3))') eig_vals(1:nconv)
           CALL putarr(fid, '/eig_vals', eig_vals(1:nconv))
           CALL putarr(fid, '/eig_vecs', eig_vecs(1:nrank,1:nconv))
           DO i=1,nconv
!!$              WRITE(*,'(/a,2(pe20.6))') '*** eigen value =', eig_vals(i)
!!$              WRITE(*,'(a/(10(1pe12.4)))') 'Real of eigen vector', &
!!$                   &                        REAL(eig_vecs(1:nrank,i))
!!$              WRITE(*,'(a/(10(1pe12.4)))') 'Imag of eigen vector', &
!!$                   &                        aimag(eig_vecs(1:nrank,i))
              WRITE(str,'(a,i3.3)') '/eig_vecs_',i
              CALL putarr(fid, TRIM(str), eig_vecs(1:nrank,i))
           END DO
           EXIT
        END IF
     END IF
  END DO  ! End of reverse commuinication loop
  IF(info .EQ. 1) THEN
     PRINT*, 'Maximum number of iterations reached!'
     PRINT*, 'IPARAM(5) =', iparam(5)
  END IF
  PRINT*, 'Number of Arnoldi iterations', iparam(3)
  tarpack = seconds() - t0
!===========================================================================
!              9.0  Epilogue
!
  WRITE(*,'(a)') '---'
  WRITE(*,'(a,1pe12.3)') 'Matrice construction time (s)  ', tmat
  IF(nlinv) THEN
     WRITE(*,'(a,1pe12.3)') 'Matrice factorization time (s) ', tfac
  END IF
  WRITE(*,'(a,1pe12.3)') 'Arpack time (s)                ', tarpack
!
  DEALLOCATE(xgrid)
  CALL destroy(mat)
  CALL destroy_sp(splx)
  CALL closef(fid)
  CALL mpi_finalize(ierr)
!
CONTAINS
  SUBROUTINE av(n,v,w)
!
    INTEGER, INTENT(in) :: n
    DOUBLE COMPLEX, INTENT(in) :: v(*)
    DOUBLE COMPLEX, INTENT(out)   :: w(*)
!
    IF(nlinv) THEN
       w(1:n) = v(1:n)
       CALL bsolve(mat,w(1:n))  ! Solve A*w = v or w=inv(A)*v
    ELSE
       w(1:n) = vmx(mat,v(1:n)) ! A*v matrix product
    END IF
  END SUBROUTINE av
!
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE PRECISION, INTENT(out) :: c(:)
!
! Weak form = Int(x*dw/dx*dt/dx)dx
!
    c(1) = x
    idt(1) = 1
    idw(1) = 1
  END SUBROUTINE coefeq
END PROGRAM main
!
!+++
SUBROUTINE meshdist(c, x, nx)
!
!   Construct an 1d  non-equidistant mesh given a
!   mesh distribution function defined in FDIST
!
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT(in) :: c(5)
  INTEGER, INTENT(iN) :: nx
  DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
  INTEGER :: nintg
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
  DOUBLE PRECISION :: a, b, dx, f0, f1, scal
  INTEGER :: i, k
!
  a=x(0)
  b=x(nx)
  nintg = 10*nx
  ALLOCATE(xint(0:nintg), fint(0:nintg))
!
!  Mesh distribution
!
  dx = (b-a)/REAL(nintg)
  xint(0) = a
  fint(0) = 0.0d0
  f1 = fdist(xint(0))
  DO i=1,nintg
     f0 = f1
     xint(i) = xint(i-1) + dx
     f1 = fdist(xint(i))
     fint(i) = fint(i-1) + 0.5*(f0+f1)
  END DO
!
!  Normalization
!
  scal = REAL(nx) / fint(nintg)
  fint(0:nintg) = fint(0:nintg) * scal
!!$  WRITE(*,'(a/(10f8.3))') 'FINT', fint
!
!  Obtain mesh point by (inverse) interpolation
!
  k = 1
  DO i=1,nintg-1
     IF( fint(i) .GE. REAL(k) ) THEN
        x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
             &   (k-fint(i))
        k = k+1
     END IF
  END DO
!
  DEALLOCATE(xint, fint)
CONTAINS
  DOUBLE PRECISION FUNCTION fdist(x)
    DOUBLE PRECISION, INTENT(in) :: x
    fdist = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
  END FUNCTION fdist
END SUBROUTINE meshdist
!+++
SUBROUTINE dismat(spl, mat)
!
!   Assembly FE matrix mat using spline spl
!
  USE bsplines
  USE mumps_bsplines
  IMPLICIT NONE
  TYPE(spline1d), INTENT(in) :: spl
  TYPE(zmumps_mat), INTENT(inout) :: mat
  INTEGER :: nrank, nx, nidbas, ngauss
  INTEGER :: i, igauss, iterm, iw, jt, irow, jcol
  DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
  DOUBLE COMPLEX                :: contrib
!
  INTEGER :: kterms                ! Number of terms in weak form
  INTEGER, ALLOCATABLE :: idert(:), iderw(:)  ! Derivative order
  DOUBLE PRECISION, ALLOCATABLE :: coefs(:)
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  CALL get_dim(spl, nrank, nx, nidbas)
  ALLOCATE(fun(0:nidbas,0:1)) ! Spline and its 1st derivative
!
!   Weak form
!
  kterms = mat%nterms
  ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms))
!
!   Gauss quadature
!
  CALL get_gauss(spl, ngauss)
  ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
  DO i=1,nx
     CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
     DO igauss=1,ngauss
        CALL basfun(xgauss(igauss), spl, fun, i)
        CALL coefeq(xgauss(igauss), idert, iderw, coefs)
        DO iterm=1,kterms
           DO jt=0,nidbas
              DO iw=0,nidbas
                 contrib = fun(jt,idert(iterm)) * coefs(iterm) * &
                      &    fun(iw,iderw(iterm)) * wgauss(igauss)
                 irow=i+iw; jcol=i+jt
                 CALL updtmat(mat, irow, jcol, contrib)
              END DO
           END DO
        END DO
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(fun)
  DEALLOCATE(xgauss, wgauss)
  DEALLOCATE(iderw, idert, coefs)
!
CONTAINS
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt))
!
! Weak form = Int(x*dw/dx*dt/dx)dx
!
    c(1) = x
    idt(1) = 1
    idw(1) = 1
  END SUBROUTINE coefeq
END SUBROUTINE dismat
!+++
SUBROUTINE disrhs(kdiff, spl, rhs)
!
!   Assenbly the RHS using spline spl
!
  USE bsplines
  IMPLICIT NONE
  INTEGER, INTENT(in) :: kdiff
  TYPE(spline1d), INTENT(in) :: spl
  DOUBLE PRECISION, INTENT(out) :: rhs(:)
  INTEGER :: nrank, nx, nidbas, ngauss
  INTEGER :: i, igauss, left
  DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
  DOUBLE PRECISION:: contrib
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  CALL get_dim(spl, nrank, nx, nidbas)
!!$  WRITE(*,'(/a, 5i3)') 'nrank, nx, nidbas =', nrank, nx, nidbas
!
  ALLOCATE(fun(nidbas+1,1)) ! needs only basis functions (no derivatives)
!
!   Gauss quadature
!
  CALL get_gauss(spl, ngauss)
!!$  WRITE(*,'(/a, i3)') 'Gauss points and weights, ngauss =', ngauss
  ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
  rhs(1:nrank) = 0.0d0
!
  DO i=1,nx
     CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
!!$     WRITE(*,'(a,i3,(8f10.4))') 'i=',i, xgauss, wgauss
     DO igauss=1,ngauss
        CALL basfun(xgauss(igauss), spl, fun, i)
        left = i-1
!!$        WRITE(*,'(a,5i4)') '>>>igauss, left', igauss, left
        contrib = wgauss(igauss)*rhseq(xgauss(igauss), kdiff)
        rhs(i:i+nidbas) = rhs(i:i+nidbas) + contrib*fun(:,1)
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(fun)
  DEALLOCATE(xgauss, wgauss)
!
CONTAINS
  DOUBLE PRECISION FUNCTION rhseq(x,k)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(in) :: k
    rhseq = k*k*x**(k-1)
  END FUNCTION rhseq
END SUBROUTINE disrhs