!>
!> @file pde1dp_cmpl_pardiso.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE pde1dp_cmpl_pardiso_mod
  USE bsplines
  USE pardiso_bsplines
  IMPLICIT NONE
!
CONTAINS
  SUBROUTINE disrhs(spl, rhs, mmode, alpha, beta)
!
  TYPE(spline1d) :: spl
  INTEGER, INTENT(in) :: mmode
  DOUBLE COMPLEX, INTENT(in) :: alpha, beta
  DOUBLE COMPLEX, INTENT(out) :: rhs(:)   
!
  INTEGER :: dim, nrank, nx, nidbas, ngauss
  INTEGER :: i, igauss, it, irow
  DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
  DOUBLE COMPLEX :: contrib
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  nrank = SIZE(rhs)
  CALL get_dim(spl, dim, nx, nidbas)
!
  ALLOCATE(fun(0:nidbas,1)) ! needs only basis functions (no derivatives)
!
!   Gauss quadature
!
  CALL get_gauss(spl, ngauss)
  ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
  rhs(:) = 0.0d0
!
  DO i=1,nx
     CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
     DO igauss=1,ngauss
        CALL basfun(xgauss(igauss), spl, fun, i)
        contrib = wgauss(igauss) * rhseq(xgauss(igauss))
        DO it=0,nidbas
           irow=MODULO(i+it-1,nx) + 1   ! Periodic BC
           rhs(irow) = rhs(irow) + contrib*fun(it,1)
        END DO
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(fun)
  DEALLOCATE(xgauss, wgauss)
CONTAINS
  DOUBLE COMPLEX FUNCTION rhseq(x)
    DOUBLE PRECISION, INTENT(in) :: x
    DOUBLE PRECISION :: arg
    arg = mmode*x
    rhseq = (mmode**2*alpha-beta)*COS(arg)
  END FUNCTION rhseq
END SUBROUTINE disrhs
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE dismat(spl, mat, alpha, beta)
!
  TYPE(spline1d) :: spl
  TYPE(zpardiso_mat) :: mat
  DOUBLE COMPLEX, INTENT(in) :: alpha, beta
!
  INTEGER :: dim, 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 COMPLEX, ALLOCATABLE :: coefs(:)
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  CALL get_dim(spl, dim, 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=MODULO(i+iw-1,nx) + 1   ! Periodic BC
                 jcol=MODULO(i+jt-1,nx) + 1
                 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 COMPLEX, INTENT(out) :: c(SIZE(idt))
!
    c(1) = alpha
    idt(1) = 1
    idw(1) = 1
!
    c(2) = -beta
    idt(2) = 0
    idw(2) = 0
  END SUBROUTINE coefeq
END SUBROUTINE dismat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  FUNCTION norm2(x)
!
!  Compute the 2-norm of complex array x
!
    IMPLICIT NONE
    DOUBLE PRECISION :: norm2
    DOUBLE COMPLEX, INTENT(in) :: x(:)
    DOUBLE PRECISION :: sum2
!
    sum2 = DOT_PRODUCT(x,x)
    norm2 = SQRT(sum2)
  END FUNCTION norm2
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE meshdist(mmode, x)
!
!   Construct a 1d  non-equidistant mesh given a
!   mesh distribution function.
!
    INTEGER, INTENT(in) :: mmode
    DOUBLE PRECISION, INTENT(inout) :: x(0:)
    INTEGER :: nx, nintg
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
    DOUBLE PRECISION :: a, b, dx, f0, f1, scal
    INTEGER :: i, k
!
    nx = SIZE(x)-1
    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
!
!  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 = 2.0 + COS(mmode*x)
    END FUNCTION fdist
  END SUBROUTINE meshdist
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
END MODULE pde1dp_cmpl_pardiso_mod
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
PROGRAM main
!
!   1D complex PDE with periodic BC
!
  USE pde1dp_cmpl_pardiso_mod
  USE futils
  USE conmat_mod
!
  IMPLICIT NONE
  TYPE(spline1d) :: splx
  TYPE(zpardiso_mat) :: mat
  TYPE(zpardiso_mat) :: newmat
  INTEGER :: kl, ku, nrank
!
  CHARACTER(len=128) :: file='pde1dp_cmpl_pardiso.h5'
  INTEGER :: fid
  INTEGER :: nx, nidbas, ngauss, mmode, npt, dim
  LOGICAL :: nlequid
  LOGICAL :: nlsym, nlherm, nlpos
  DOUBLE PRECISION, PARAMETER :: pi=3.1415926535897931d0
  DOUBLE PRECISION :: dx
  DOUBLE COMPLEX :: alpha, beta
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, x, err
  DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE :: sol, rhs, bcoef
  DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE :: newsol, arow
  DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE :: solcal, solana
  DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
  DOUBLE PRECISION :: err_norm
  INTEGER :: i
!
  NAMELIST /newrun/ nx, nidbas, ngauss, nlequid, alpha, beta, mmode, npt, &
       &            nlsym, nlherm, nlpos
!===========================================================================
!              1.0 Prologue
!
!   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
  nlequid = .TRUE.   ! Use exact sol. as mesh dist. function if .FALSE.
  mmode = 1          ! Fourier mode
  alpha = (1.0, 1.0) ! Complex "diffusion"
  beta  = 1.0
  npt = 100
  nlsym = .TRUE.    ! Is matrice symmetric
  nlherm = .FALSE.  ! Is matrice hermitian
  nlpos = .TRUE.    ! and positive definite ?
!
  READ(*,newrun)
  WRITE(*,newrun)
!
!   Define grid on x axis
!
  ALLOCATE(xgrid(0:nx))
  dx = 2.d0*pi/REAL(nx,8)
  xgrid = (/ (i*dx,i=0,nx) /)
  IF( .NOT. nlequid ) THEN
     CALL meshdist(mmode, xgrid)
  END IF
  WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
!
!   Create hdf5 file
!
  CALL creatf(file, fid, 'PDE1DP Result File', real_prec='d')
  CALL attach(fid, '/', 'NX', nx)
  CALL attach(fid, '/', 'NIDBAS', nidbas)
  CALL putarr(fid, '/xgrid', xgrid)
!===========================================================================
!              2.0 Discretize the PDE
!
!   Set up periodic spline
!
  CALL set_spline(nidbas, ngauss, xgrid, splx, period=.TRUE.)
  WRITE(*,'(a,l6)') 'nlequid =', nlequid
  nrank = nx       ! Rank of the FE matrix
!
!   FE matrix assembly
!
  CALL init(nrank, 2, mat, nlsym=nlsym, nlherm=nlherm, nlpos=nlpos)
  CALL get_dim(splx, dim)
  WRITE(*,'(/a,4i6)') 'nrank, dim', nrank, dim
!!$  CALL dismat(splx, mat, alpha, beta)
  CALL conmat(splx, mat, coefeq)
!
!   RHS assembly
!  
  ALLOCATE(rhs(nrank), sol(nrank))
  CALL disrhs(splx, rhs, mmode, alpha, beta)
!
  CALL putarr(fid, '/rhs', rhs)
!
!   Factor and solve
!
  WRITE(*,'(a/(10i6))') 'iparm', mat%p%iparm
  CALL factor(mat)
  CALL putmat(fid,'/MAT', mat)
  CALL bsolve(mat, rhs, sol)
  CALL putarr(fid, '/sol', sol)
  WRITE(*,'(/a,i8)') 'Number of non-zeros of matrix       = ',get_count(mat)
  WRITE(*,'(a,i8)')  'Number of nonzeros in factors of A  = ',mat%p%iparm(18)
!
!   Compute residue
!
  WRITE(*,'(a,1pe12.4)') 'Residue =', norm2(vmx(mat,sol)-rhs)
!===========================================================================
!              3.0   Check solution
!
!   Exact solution
  ALLOCATE(x(0:npt), solcal(0:npt), solana(0:npt), err(0:npt))
  dx=2.0d0*pi/REAL(npt,8)
  x = (/ (i*dx, i=0,npt) /)
  solana = COS(mmode*x)
!
!   Prolongate solution using periodicity
!
  ALLOCATE(bcoef(dim))
  bcoef(1:nrank) = sol(1:nrank)
  DO i=nrank+1,dim
     bcoef(i) = bcoef(MODULO(i-1,nrank)+1)
  END DO
!
!   Interpolate field
!
  CALL gridval(splx, x, solcal, 0, bcoef)
!
  err = ABS(solcal-solana)
  CALL putarr(fid, '/x', x)
  CALL putarr(fid, '/solana', solana)
  CALL putarr(fid, '/solcal', solcal)
  CALL putarr(fid, '/err', err)
!
!   Compute discretization error norm by Gauss integration
!
  err_norm=0.0
  ALLOCATE(xgauss(ngauss), wgauss(ngauss))
  DO i=1,nx
     CALL get_gauss(splx, ngauss, i, xgauss, wgauss)
     CALL gridval(splx, xgauss, solcal(1:ngauss), 0)
     solana(1:ngauss) = COS(mmode*xgauss)
     err(1:ngauss) = DOT_PRODUCT(solana(1:ngauss)-solcal(1:ngauss), &
          &          solana(1:ngauss)-solcal(1:ngauss))
     err_norm = err_norm + SUM(wgauss*err(1:ngauss))
  END DO
  err_norm = SQRT(err_norm)
  WRITE(*,'(a,1pe12.3)') 'Discretization error ',  err_norm
!
!===========================================================================
!              4.0 Test of getrow/putrow, getcol/putcol and mcopy
!
  CALL init(nrank, 2, newmat, nlsym=nlsym, nlherm=nlherm, nlpos=nlpos)
  ALLOCATE(arow(nrank), newsol(nrank))
!
  DO i=1,nrank
     CALL getrow(mat, i, arow)
     CALL putrow(newmat, i, arow)
  END DO
  CALL factor(newmat)
  CALL bsolve(newmat, rhs, newsol)
  WRITE(*,'(/a)') 'putrow/getrow ...'
  WRITE(*,'(a,1pe12.4)') 'Residue =', norm2(vmx(newmat,newsol)-rhs)
  WRITE(*,'(a,1pe12.3)') 'Error ',  norm2(sol-newsol)
!
  DO i=1,nrank
     CALL getcol(mat, i, arow)
     CALL putcol(newmat, i, arow)
  END DO
  CALL factor(newmat)
  CALL bsolve(newmat, rhs, newsol)
  WRITE(*,'(/a)') 'putcol/getcol ...'
  WRITE(*,'(a,1pe12.4)') 'Residue =', norm2(vmx(newmat,newsol)-rhs)
  WRITE(*,'(a,1pe12.3)') 'Error ',  norm2(sol-newsol)
!
  CALL clear_mat(newmat)
  CALL mcopy(mat, newmat)
  WRITE(*,'(/a)') 'mcopy ...'
  newmat%val = (1000.0d0,0.0d0)*newmat%val
  CALL factor(newmat)
  CALL bsolve(newmat, rhs, newsol)
  WRITE(*,'(a)') 'Backsolve the new system'
  WRITE(*,'(a,1pe12.4)') 'Norm =', norm2(newsol)
!
  WRITE(*,'(a)') 'Destroy NEWMAT ...'
  CALL destroy(newmat)
!
  CALL bsolve(mat, rhs, sol)
  WRITE(*,'(/a)') 'Backsolve the old system'
  WRITE(*,'(a,1pe12.4)') 'Norm =', norm2(sol)
!
  WRITE(*,'(a)') 'Destroy MAT ...'
  CALL destroy(mat)
!!$!
!!$  WRITE(*,'(/a)') 'Should crash since NEWMAT is gone!'
!!$  CALL bsolve(newmat, rhs, newsol)
!!$  WRITE(*,'(a)') 'Backsolve the new system'
!!$  WRITE(*,'(a,1pe12.4)') 'Norm =', norm2(newsol)
!===========================================================================
!              9.0  Clean up
!
  DEALLOCATE(x, solcal, solana, err)
  DEALLOCATE(xgauss, wgauss)
  DEALLOCATE(bcoef)
  DEALLOCATE(xgrid)
  DEALLOCATE(rhs, sol)
  DEALLOCATE(arow, newsol)
  CALL destroy_sp(splx)
  CALL closef(fid)  
CONTAINS
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE COMPLEX, INTENT(out) :: c(:)
!
    c(1) = alpha
    idt(1) = 1
    idw(1) = 1
!
    c(2) = -beta
    idt(2) = 0
    idw(2) = 0
  END SUBROUTINE coefeq
END PROGRAM main